Flows of multiple fluid phases are common in many subsurface reservoirs. Numerical simulation of these flows can be
challenging and computationally expensive. Dynamic adaptive mesh optimisation and related approaches, such as adaptive
grid refinement can increase solution accuracy at reduced computational cost. However, in models or parts of the model
domain, where the local Courant number is large, the solution may propagate beyond the region in which the mesh is
refined, resulting in reduced solution accuracy, which can never be recovered. A methodology is presented here to modify
the mesh within the non-linear solver. The method allows efficient application of dynamic mesh adaptivity techniques even
with high Courant numbers. These high Courant numbers may not be desired but a consequence of the heterogeneity of the
domain. Therefore, the method presented can be considered as a more robust and accurate version of the standard dynamic
mesh adaptivity techniques.