Numerical models of magmatic hydrothermal systems have become powerful tools for linking surface and seafloor observations to chemical and fluid-dynamic processes at depth. This task requires resolving multi-phase flow over large distances of several kilometers, a wide range of pressure (p) and temperature (T) conditions, and over timescales of several thousands of years. The key numerical challenge is that realistic simulations have to consider the high nonlinearity and strong coupling of the governing conservation equations for mass and energy, while also being numerically efficient so that the required spatial and temporal scales can be resolved. Here we outline possible solutions to this problem by evaluating different implementation strategies and presenting a numerical scheme for fully coupled accurate and efficient flow solutions. The general scheme, based on the Newton-Raphson (NR) method, is presented for the simplified case of 2-D pure water convection and uses a control volume discretization on unstructured meshes. We find that the presented techniques significantly reduce the computational effort with respect to sequential/decoupled schemes. Key to this is a theta-time-differencing method for better accuracy, stability and convergence behavior of the NR-iterations, as well as improvements regarding upwinding. These features make the presented methods useful for coupled simulations of magmatic hydrothermal systems and a potential basis for future 3-D multi-phase codes.