We study the evolution of hypersurfaces in spacetime initial data sets by their null mean curvature. A theory of weak solutions is developed using the level-set approach. Starting from an arbitrary mean convex, outer untapped hypersurface $\partial\Omega_0$, we show that there exists a weak solution to the null mean curvature flow, given as a limit of approximate solutions that are defined using the $\varepsilon$-regularization method. We show that the approximate solutions blow up on the outermost MOTS and the weak solution converges (as boundaries of finite perimeter sets) to a generalized MOTS.
