Coastal zones are vulnerable to both erosion and flood risk, which can be assessed using coupled hydro-
morphodynamic models. However, the use of such models as decision support tools suffers from a high degree
of uncertainty, due to both incomplete knowledge and natural variability in the system. In this work, we show
for the first time how the multilevel Monte Carlo method (MLMC) can be applied in hydro-morphodynamic
coastal ocean modelling, here using the popular model XBeach, to quantify uncertainty by computing statistics
of key 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. Several theoretical and real-world coastal
zone case studies are considered here, for which output variables that are key to the assessment of flood and
erosion risk, such as wave run-up height and total eroded volume, are estimated. We show that MLMC can
significantly reduce computational cost, resulting in speed up factors of 40 or greater compared to a standard
Monte Carlo approach, whilst keeping the same level of accuracy. Furthermore, a sophisticated ensemble
generating technique is used to estimate the cumulative distribution of output variables from the MLMC output.
This allows for the probability of a variable exceeding a certain value to be estimated, such as the probability
of a wave run-up height exceeding the height of a seawall. This is a valuable capability that can be used to
inform decision-making under uncertainty