  Mathematical models for flow through porous media typically enjoy the
so-called maximum principles, which place bounds on the pressure field. It is
highly desirable to preserve these bounds on the pressure field in predictive
numerical simulations, that is, one needs to satisfy discrete maximum
principles (DMP). Unfortunately, many of the existing formulations for flow
through porous media models do not satisfy DMP. This paper presents a robust,
scalable numerical formulation based on variational inequalities (VI), to model
non-linear flows through heterogeneous, anisotropic porous media without
violating DMP. VI is an optimization technique that places bounds on the
numerical solutions of partial differential equations. To crystallize the
ideas, a modification to Darcy equations by taking into account
pressure-dependent viscosity will be discretized using the lowest-order
Raviart-Thomas (RT0) and Variational Multi-scale (VMS) finite element
formulations. It will be shown that these formulations violate DMP, and, in
fact, these violations increase with an increase in anisotropy. It will be
shown that the proposed VI-based formulation provides a viable route to enforce
DMP. Moreover, it will be shown that the proposed formulation is scalable, and
can work with any numerical discretization and weak form. Parallel scalability
on modern computational platforms will be illustrated through strong-scaling
studies, which will prove the efficiency of the proposed formulation in a
parallel setting. Algorithmic scalability as the problem size is scaled up will
be demonstrated through novel static-scaling studies. The performed
static-scaling studies can serve as a guide for users to be able to select an
appropriate discretization for a given problem size.
