  In this paper, we construct global-in-time forward and backward Lagrangian
flow maps along the pressure gradient generated by weak solutions of the Porous
Media Equation. The main difficulty is that when the initial data has compact
support, it is well-known that the pressure gradient is not a BV function.
Thus, the theory of regular Lagrangian flows cannot be applied to construct the
flow maps. To overcome this difficulty, we develop a new argument that combines
Aronson-B\'enilan type estimates with the quantitative Lagrangian flow theory
of Crippa and De Lellis to show that certain doubly logarithmic quantities
measuring the stability of flow maps do not blow up fast enough to prevent
compactness. Our arguments are sufficiently flexible to handle the Hele-Shaw
limit and a multispecies generalization of the Porous Media Equation where the
equation is replaced by a coupled hyperbolic-parabolic system of reaction
diffusion equations. As one application of our flow maps, we are able to
construct solutions where different species cannot mix together if they were
separated at initial time.
