Dense granular flows can spontaneously self-channelise by forming a pair of parallelsided static levees on either side of a central flowing channel. This process prevents lateral spreading and maintains the flow thickness, and hence mobility, enabling the grains to run-out considerably further than a spreading flow on shallow slopes. Since levees commonly form in hazardous geophysical mass flows, such as snow avalanches, debrisflows, lahars and pyroclastic flows, this has important implications for risk management in mountainous and volcanic regions. In this paper an avalanche model that incorporates frictional hysteresis, as well as depth-averaged viscous terms derived from the μ(I)- rheology, is used to quantitatively model self-channelisation and levee formation. The viscous terms are crucial for determining a smoothly varying steady-state velocity profile across the flowing channel, which has the important property that it does not exert any shear-stresses at the levee-channel interfaces. For a fixed mass flux, the resulting boundary value problem for the velocity profile also uniquely determines the width and height of the channel, and the predictions are in very good agreement with existing experimental data for both spherical and angular particles. It is also shown that in the absence of viscous (second-order gradient) terms, the problem degenerates, to produce plug-flow in the channel with two frictionless contact discontinuities at the levee-channel margins. Such solutions are not observed in experiments. Moreover, the steady-state inviscid problem lacks a thickness or width selection mechanism and consequently there is no unique solution. The viscous theory is therefore a significant step forward. Fully time-dependent numerical simulations to the viscous model are able to quantitatively capture the process in which the flow self-channelizes and show how the levees are initially emplaced behind the flow head. Both experiments and numerical simulations show that the height and width of the channel are not necessarily fixed by these initial values, but respond to changes in the supplied mass flux, allowing narrowing and widening of the channel long after the initial front has passed by. In addition, below a critical mass flux the steadystate solutions become unstable and time-dependent numerical simulations are able to capture the transition to periodic erosion-deposition waves observed in experiments.