“…A number of different numerical methods have been developed to handle the gaseousaqueous phase transitions in multi-phase multi-component porous media models, e.g., primary variable switching (PVS) schemes [47,8], method of negative saturations [41], method of persistent variables [40,26], and non-linear complementary constraints (NCP) approaches P r e p r i n t An all-at-once Newton strategy for methane hydrate reservoir models 3 [36,34,3,5]. In the most widely used gas hydrate reservoir simulators, e.g., TOUGH-Hydrate [38], HYRES-C [29,28], STOMP-HYD [45], HRS [39], etc., the gaseous-aqueous phase transitions are handled using the PVS schemes, where the choice of the primary variables is adapted locally to the phase state. However, due to the strong coupling and nonlinearities, the phase states in gas hydrate models tend to switch back and forth rapidly, and this often leads to spurious oscillation and a drastic reduction in time step size, in the extreme case, even to a breakdown of the numerical algorithm.…”