The potential benefits and risks of natural methane hydrate occurrences have raised the need to understand the processes governing hydrate formation and dissociation over the last decades, and models of increasing complexity have been developed for this purpose. We propose a formulation of a multi-dimensional methane hydrate model that couples the established chemistry of hydrate formation to the fluid flow in marine sediments undergoing compaction. The numerical model applies the Finite Volume Method to construct a segregated solver for the coupled system. The solution of the sequence of individual processes is stabilised with an adaptive Picard iteration in each time step. We implement the model based on the OpenFOAM library and extend the functionality by generalising the diffusion term to take into account variable porosity of the sediment matrix. The advantage of this robust and efficient scheme is that the formulation is conservative by construction, allowing an accurate solution of mass transport when strong concentration gradients develop in the vicinity of phase boundaries. We validate the model using data from the Blake Ridge hydrate province. Good agreement is found for the pore water chlorinity, a proxy for hydrate formation, except for regions where heterogeneous hydrate formation results in highly variable measurements.