The flow of fluids within porous rocks is an important process with numerous applications in Earth sciences. Modeling the compaction-driven fluid flow requires the solution of coupled nonlinear partial differential equations that account for the fluid flow and the solid deformation within the porous medium. Despite the nonlinear relation of porosity and permeability that is commonly encountered, natural data show evidence of channelized fluid flow in rocks that have an overall layered structure. We consider a system of nonlinear PDEs for porosity and effective pressure, based on a poroviscoelastic model, which describes such phenomena. We first discuss well-posedness of this PDE problem, which had been established in the literature only for initial porosities of high Sobolev smoothness. We present several results for porosities of low regularity, including cases with jump discontinuities that are of particular interest in geological applications since layers of different rock types routinely have discontinuous hydraulic and mechanical properties. Then we turn to results on an adaptive numerical method, which is based on a fixed-point scheme inspired by the analysis, combined with a space-time least-squares formulation. This yields an appropriate treatment of discontinuities and enables spatially varying time steps, which are required for efficient approximations of the strongly spatially and temporally localized features of solutions. We numerically show its quasi-optimality even for discontinuous porosity distributions. Furthermore, we show numerical results on a preconditioner for the space-time least-squares method to allow for efficient iterative solvers also in higher dimensions. Lastly, our approach enables a straightforward coupling to models of mass transport for chemical trace elements as the entire evolution history is stored efficiently. Our results show the influence of different kinds of layering in the development of fluid-rich channels and mass transport.