The risk from erosion and flooding in the coastal zone has the potential to increase in a changing climate. The development and use of coupled hydro-morphodynamic models is therefore becoming an ever higher priority. However, their use as decision support tools suffers from the high degree of uncertainty associated with them, due to incomplete knowledge as well as natural variability in the system. Here we show for the first time how the multilevel Monte Carlo method (MLMC) can be applied to hydro-morphodynamic models, in this case XBeach, to quantify uncertainty by computing statistics of output variables given uncertain input parameters. MLMC accelerates the Monte Carlo approach through the use of a hierarchy of models with different levels of resolution. A variety of theoretical and real-world coastal zone case studies are considered, for which output variables that are important to the assessment of flood and erosion risk are estimated, such as wave run-up height and total eroded volume. We show that MLMC can significantly reduce computational cost, resulting in speed up factors of 40 or greater compared to a simple Monte Carlo approach, whilst still maintaining the same level of accuracy. Furthermore, MLMC is used to estimate the cumulative distribution of these output variables for given uncertain parameters. This allows the risk of a variable exceeding a certain value to be calculated, for example the risk of the wave run-up height exceeding the height of a physical structure such as a seawall; this is a useful capability to inform decision-making processes.