We are solving a problem of salinisation of coastal aquifers. As a test case example, we consider the Henry saltwater intrusion problem. Since porosity, permeability and recharge are unknown or only known at a few points, we model them using random fields and random variables. The Henry problem describes a two-phase flow and is non-linear and time-dependent. The solution to be found is the expectation of the salt mass fraction, which is uncertain and time-dependent. To estimate this expectation, we use the well-known multilevel Monte Carlo (MLMC) method. The MLMC method takes just a few samples on computationally expensive (fine) meshes and more samples on cheap (coarse) meshes. Then, by building a telescoping sum, the MLMC method estimates the expected value at a much lower computational cost than the classical Monte Carlo method. The deterministic solver used here is the well-known parallel and scalable ug4 solver.