Many problems in biology involve gels which are mixtures composed of a polymer network permeated by a ﬂuid solvent (water). The two-ﬂuid model is a widely used approach to described gel mechanics, in which both network and solvent coexist ateachpoint of space andtheir relativeabundanceis describedby their volume fractions. Each phase is modeled as a continuum with its own velocity and constitutive law. In some biological applications, free boundaries separate regions of gel and regions of pure solvent, resulting in a degenerate network momentum equation where the network volume fraction vanishes. To overcome this difﬁculty, we develop a regularization method to solve the two-phase gel equations when the volume fraction of one phase goes to zero in part of the computational domain. A small and constant network volume fraction is temporarily added throughout the domain in setting up the discrete linear equations and the same set of equation is solved everywhere. These equations are very poorly conditioned for small values of the regularization parameter, but the multigrid-preconditioned GMRES method we use to solve them is efﬁcient and produces an accurate solution of these equations for the full range of relevant regularization parameter values.