This model of the subglacial drainage system simulates the pressurised flow of water at the ice-bed interface of glaciers and ice sheets. It includes both distributed and channelized water flow. Notably the model determines the geometry of the channel network as part of the solution. The resulting channel network is similar to subaerial stream networks with channels carving out hydraulic potential "valleys". However, there are some pronounced differences to subaerial drainage, for example that the time for a network to form (and decay) is on the order of weeks to months; or that, channels originating at point sources can lie on ridges of the hydraulic potential. The model employs a novel finite element approach to solve the parabolic equations for the hydraulic potential simultaneously on the 1D channel network and 2D distributed system.