Marine gas hydrate systems are characterized by highly dynamic transport-reaction processes in an essentially water-saturated porous medium that are coupled to thermodynamic phase transitions between solid gas hydrates, free gas and dissolved methane in the aqueous phase. These phase transitions are highly nonlinear and strongly coupled, and cause the mathematical model to rapidly switch the phase states and pose serious convergence issues for the classical Newton's method. One of the common methods of dealing with such phase transitions is the primary variable switching (PVS) method where the choice of the primary variables is adapted locally to the phase state 'outside' the Newton loop. In order to ensure that the phase states are determined accurately, the PVS strategy requires an additional iterative loop, which can get quite expensive for highly nonlinear problems. For methane hydrate reservoir models, the PVS method shows poor convergence behaviour and often leads to extremely small time step sizes. In order to overcome this issue, we have developed a nonlinear complementary constraints method (NCP) for handling phase transitions, and implemented it within a non-smooth Newtons linearization scheme using an active-set strategy. Here, we present our numerical scheme and show its robustness through field scale applications based on the highly dynamic geological setting of the Black Sea.