In this paper, we describe a numerical algorithm for the self-consistent simulations of surface water and sediment dynamics. The method is based on the original Lagrangian-Eulerian CSPH-TVD approach for solving the Saint-Venant and Exner equations, taking into account the physical factors essential for the understanding of the shallow water and surface soil layer motions, including complex terrain structure and its evolution due to sediment transport. Additional Exner equation for sediment transport has been used for the numerical CSPH-TVD scheme stability criteria definition. By using OpenMP-CUDA and GPUDirect technologies for hybrid computing systems (supercomputers) with several graphic coprocessors (GPUs) interacting with each other via the PCI-E / NVLINK interface we also develop a parallel numerical algorithm for the CSPH-TVD method. The developed parallel version of the algorithm demonstrates high efficiency for various configurations of Nvidia Tesla CPU + GPU computing systems. In particular, maximal speed up is 1800 for a system with four C2070 GPUs compare to the serial version for the CPU. The calculation time on the GPU V100~(Volta architecture) is reduced by 95 times compared to the GPU C2070~(Fermi architecture).