Many open problems in the Earth sciences can only be understood by modelling the porous flow of melt through a viscously deforming solid rock matrix. However, the system of equations describing this process becomes mathematically degenerate in the limit of vanishing melt fraction... Numerical methods that do not consider this degeneracy or avoid it solely by regularising specific material properties generally become computationally expensive as soon as the melt fraction approaches zero in some part of the domain. Here, we present a new formulation of the equations for coupled magma/mantle dynamics that addresses this problem, and allows it to accurately compute large-scale 3-D magma/mantle dynamics simulations with extensive regions of zero melt fraction. We achieve this by rescaling one of the solution variables, the compaction pressure, which ensures that for vanishing melt fraction, the equation causing the degeneracy becomes an identity and the other two equations revert to the Stokes system. This allows us to split the domain into two parts, and to solve the coupled magma/mantle dynamics system only in cells where melt is present. We have implemented this formulation in the open source geodynamic modelling code ASPECT and illustrate the improved performance compared to the previous three-field formulation, showing numerically that the new formulation is optimal in terms of problem size and only minimally sensitive to model parameters. Beyond that, we demonstrate the applicability to realistic problems by showing large-scale 2-D and 3-D models of mid-ocean ridges with complex rheology. Hence, we believe that our new formulation and its implementation in ASPECT will prove a valuable tool for studying the interaction of melt segregating through and interacting with a solid host rock in the Earth and other planetary bodies using high-resolution, three-dimensional simulations. (read more)