Two fluid flow in porous medium systems is an important application in many different areas of science and engineering.  Overwhelmingly, it is necessary to mathematically model the behavior of applications of concern at an averaged scale where the juxtaposed position of the phases is not resolved in detail.  This length scale is called the macroscale and the traditional model that is used nearly universally was formulated phenomenologically nearly 100 years ago. Since that time, considerable, and important, work has been done to close this model, advance more efficient numerical methods to solve the resultant equations, and analyze mathematical aspects of the model behavior.   This considerable work notwithstanding, many important open issues remain with this traditional model, including hysteretic behavior of the closure relations, lack of connection between the microscale and the macroscale, absence of explicit dependence upon variables well known to be important (contact angles, interfacial tensions, curvatures, etc), and the absence of thermodynamic constraints.  We report on a sustained effort to resolve these theoretical shortcomings and to formulate a new generation of models for this important class of applications.  A summary of the theoretical approach based upon the thermodynamically constrained averaging theory (TCAT) is discussed, a hierarchy of models is formulated, and an example model instance is examined in detail.   The problem of model closure and validation is considered.  Relying upon notions from integral geometry, we formulate a hysteretic-free state equation that applies under both equilibrium and dynamic conditions. We show that this equation is essentially exact by comparing to high-resolution simulations for a wide variety of systems.  We extend the notion of a state equation to resistance coefficients, and we show promising results for the removal of hysteresis from common relative permeability relations as well.  We summarize recent results to derive evolution equations for curvatures, and we assemble the various components to reveal a complete, closed model.