A mathematical formulation is presented for 2-D depth average unsteady flow, sediment transport, and bed evolution computation in natural watercourses. The model is developed in a general curvilinear coordinate system to characterize the longitudinal and traverse movement; it is designed to study small reservoirs, rivers reaches and estuaries, wherever depth averaging is appropriate. The model is built in a modular form, therefore different flow resistance and sediment transport criteria can be chosen. Sediment mixtures in natural watercourses are represented through a suitable number of discrete size classes, of which some of them may be subject to bedload transport depending on prevailing local hydraulic conditions. The system of sediment equations is described for: n-size class mass-conservation equations, describing bedload transport and size-fraction change, and a global mass-conservation equation for all bed material, describing the change in bed surface level. The system of equations at any computational point is solved by an implicit finite difference scheme. Starting from a computed free-surface flow, the sediment equations system is computed and coupled with the water flow in an iterative way in order to take into account the interaction between the flow field and changes in both, bed elevation and bed surface size distribution. The mathematical formulation is with experimental data from measurements of the hydraulics laboratory of the Engineering Institute.