The mathematical model of shear shallow water flows of constant density is studied. This is a 2D hyperbolic non-conservative system of equations that is mathematically equivalent to the Reynolds-averaged model of barotropic turbulent flows. The model has three families of characteristics corresponding to the propagation of surface waves, shear waves and average flow (contact characteristics). The system is non-conservative: for six unknowns (the fluid depth, two components of the depth averaged horizontal velocity, and three independent components of the symmetric Reynolds stress tensor) one has only five conservation laws (conservation of mass, momentum, energy and mathematical ‘entropy’). A splitting procedure for solving such a system is proposed allowing us to define a weak solution. Each split subsystem contains only one family of waves (either surface or shear waves) and contact characteristics. The accuracy of such an approach is tested on 2D analytical solutions describing the flow with linear with respect to the space variables velocity, and on the solutions describing 1D roll waves. The capacity of the model to describe the full transition scenario as commonly seen in the formation of roll waves: from uniform flow to 1D roll waves, and, finally, to 2D transverse ‘fingering’ of the wave profiles, is shown.