Shallow-water type models are commonly used in tsunami simulations. These models contain uncertain parameters like the ratio of densities of layers, friction coefficient, fault deformation, etc. These parameters are modeled statistically and quantifying the resulting solution uncertainty (UQ) is a crucial task in geophysics. We propose a paradigm for UQ that combines the recently developed path-conservative spatial discretizations efficiently implemented on cluster of GPUs, with the recently developed Multi-Level Monte Carlo (MLMC) statistical sampling method and provides a fast, accurate and computationally efficient framework to compute statistical quantities of interest. Numerical experiments, including realistic simulations in real bathymetries, are presented to illustrate the robustness of the proposed UQ algorithm.