Darcy, Forchheimer, Brinkman and Richards: classical hydromechanical equations and their significance in the light of the TPM

In hydromechanical applications, Darcy, Brinkman, Forchheimer and Richards equations play a central role when porous media flow under saturated and unsaturated conditions has to be investigated. While Darcy, Brinkman, Forchheimer and Richards found their equations mainly on the basis of flow observations in field and laboratory experiments, the modern Theory of Porous Media allows for a scientific view at these equations on the basis of precise continuum mechanical and thermodynamical investigations. The present article aims at commenting the classical equations and at deriving their counterparts by the use of the thermodynamical consistent Theory of Porous Media. This procedure will prove that the classical equations are valid under certain restrictions and that extended equations exist valid for arbitrary cases in their field.


Introduction
The present article aims at discussing the validity of classical hydromechanical equations describing porous media flow for saturated and unsaturated porous materials. Thus, the article is somehow in between a review article and an article presenting original research. As far as the author is aware, the famous Darcy law is widely used as a constitutive assumption for the description of mostly fully saturated pore flow conditions in porous media, rather than as the result of a continuum mechanical and thermodynamical investigation of a biphasic material of solid and fluid. In the same way, the Brinkman or Darcy-Brinkman equation as well as the Forchheimer or Darcy-Forchheimer equation is handled as extensions of Darcy's law either by the introduction of frictional forces or by including inertia-or tortuosity-based terms added to the drag force. In contrast to these two equations, the Richards equation has been formulated for fluid flow in unsaturated media by combining the Darcy equation with the continuity equation of the pore fluid.
In the following chapters, the classical approaches will be addressed, before biphasic and triphasic media are considered with the goal to find general relations that can be reduced to the classical equations. By this approach, it will easily be seen, what assumptions have to be made, such that the classical equations become evident. As the classical equations do not include thermal effects, the following investigations will be carried out under the assumption of isothermal conditions.

Darcy 1856
In 1856, Henry Darcy [1] published an extended treatise on the public water supply for the city of Dijon, compare Fig. 1, where he also included the results of various experiments that guided him to formulate his famous filter law yielding "Il paraît donc que, pour un sable de même nature, on peut admettre que le volume débité est proportionnel à la charge et en raison inverse de l'épaisseur de la couche traversé". This led to the equation where grad (·) = ∂ (·)/∂x with x as the location vector. Darcy's empirical equation given here in its basic one-dimensional (1-d) and in a modern three-dimensional (3-d) notation expresses the filter velocity v F V of the pore fluid as a linear function of the gradient of the hydraulic head h. While v F V is given in m/s, h is measured in m and the hydraulic conductivity k F in m 3 /(m 2 s), often reduced to m/s, respectively. Note that the hydraulic head is usually given as cf. Figure 2, where γ F R = ρ F R g is the intrinsic or realistic specific weight of the fluid, while ρ F R is the intrinsic fluid density that is constant as a result of the assumption of fluid incompressibility. Furthermore, p F R is the intrinsic fluid pressure exceeding the atmospheric pressure. Finally, g = |g| is the norm of the gravitation vector g = −g e 3 . Note in passing that h can also be understood as the so-called water potential splitting in a pressure-driven and a gravitation-driven term. Dividing the gravitation potential U g by the specific weight yields with U 0 = 0 Therein, p F R /γ F R is the pressure head, sometimes also called the matrix head, while x 3 is the geodetic head.
As an alternative to (1), Darcy's filter law can also be formulated by the use of the intrinsic permeability K S given in m 2 , such that While K S only depends on the permeability properties of the porous medium, k F additionally depends on the specific weight γ F R and the intrinsic shear viscosity μ F R of the fluid. Introducing a specific permeability k through with the dimension m 3 s/kg, Darcy's law yields v F V = −k (grad p F R + ρ F R g).
2. 2 Forchheimer 1901 Although Darcy's filter law was accepted as a more or less general equation describing the movement of water through soil, there have been several issues requiring extensions. Especially, it has been found that the proportionality between the gradient of the hydraulic head and the filter velocity of the pore fluid did not persist for higher flow rates. This phenomenon has been the subject of many experimental investigations and theoretical speculations, where microinertia terms or tortuosity properties came into play. In the original article, Forchheimer [2] argues on the basis of experimental findings on coarse sands that the hydraulic conductivity shrinks when the hydraulic gradient increases, such that the proportionality between v F V and grad h does not hold. Based on this, Forchheimer suggested to modify Darcy's law (1) 1 by the addition of a quadratic velocity term via with a = 1/k F and b as a non-Darcian flow coefficient. Apart of the parabola (7), Forchheimer additionally suggested a cubic version of (7) and a potential function −∂h/∂ x = m (v F V ) n . Generalising the one-dimensional equation (7) towards three dimensions yields In the younger literature, compare, for example, Markert [3], b is set to b = 1/b S , where b S with the dimension (m/s) 2 is a sort of solid tortuosity parameter including the inherent irregularities of the pore channels in the description.
In extension of Forchheimer's work, (8) can be solved for the filter velocity where |v F V | still needs to be expressed by |grad h|. This can easily been done by taking the norm of both sides of (8): As a quadratic equation in |v F V |, the solution of (10) is Based on the fact that |v F V | has to vanish when |grad p F R | vanishes, only |v F V | 1 is valid, such that (9) can be solved as a function of grad h alone: In those cases, where the hydraulic gradient should be expressed by the pressure gradient, grad h has to be substituted by where (2) has been used. With the aid of this equation, (12) where k F has been substituted by K S trough (5) and b S has been substituted by the tortuosity parameter B S with dimension m through b S = B S |g|, thus combining the tortuosity effect with the norm of the gravitation g.
In case of negligible body forces, (14) reduces to Note that in equations (12), (14) and (15), the terms K F and K F F can be understood as hydraulic conductivities that, in contrast to k F , also consider the relative importance of tortuosity-induced inertia forces, cf. Markert [3]. Furthermore, note that using the intrinsic measures K S and B S is generally preferable, as they are independent of the properties of the permeating fluid.

Brinkman 1949
As Forchheimer, Brinkman [4] was looking for a modification of Darcy's law for the description of, as the title of [4] says, "A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles". Starting with Darcy's law (4), where the gravitational forces have been neglected, such that respectively, Brinkman considered this equation as not appropriate to solve the problem as there is no definition of a viscous force. Following this, he looked at the Navier-Stokes equations of a viscous, incompressible fluid yielding where v F is the fluid velocity or the seepage velocity of a pore fluid permeating a non-deforming skeleton, respectively. The seepage velocity is related to the filter velocity through v F V = n F v F . Therein, n F is the fluid volume fraction or, in case of a fully saturated porous solid, the porosity. Furthermore, div (·) is the divergence operator corresponding to grad (·), while Δ (·) = div grad (·) is the Laplace operator. Note that (17) 1 is the continuity equation of a free fluid flow, while (17) 2 is the momentum balance of a viscous fluid, where inertia and gravitational forces have been neglected. Stating that the Navier-Stokes equations combine partial differential equations of first and second order "making it impossible to formulate rational boundary equations", Brinkman suggested to combine the Darcy and the Navier-Stokes equations via whereμ F R was assumed to be possibly different from μ F R . As far as the author is aware, Brinkman did not recognise the difference between v F and v F V . Writing the above equation as Brinkman did in terms of v F , Brinkman would have found that where μ F = n F μ F R is the partial shear viscosity of the pore fluid, and Brinkman's assumption that there is a difference between the two viscosity parameters is obvious. As Brinkman stated in his article, it is furthermore seen from (19) that "this equation has the advantage of approximating (16) 2 for low values of K S and (17) 2 for high values of K S ". Concerning Brinkman's equation and the consequences of its individual terms, also compare the work by Auriault and colleagues [5,6], where the second term in (19) was called a corrector to the original Darcy equation.

Richards 1931
In contrast to Darcy, Forchheimer and Brinkman, who considered a fluid-saturated porous medium, Richards [7] was interested in a partially saturated material. Based on earlier work by Buckingham [8] and Richardson [9], he combined Darcy's law with the continuity equation of the pore fluid. In a partially saturated material, the continuity equation of the pore liquid is obtained from the mass balance of the pore liquid yielding in an Eulerian setting under the assumption of liquid incompressibility Therein, n L is the volumetric fraction of the pore liquid compared to the overall volume of the porous medium of solid, liquid and gas. Taking the divergence of Darcy's law in the version of (1) 2 results in where Richards substituted the hydraulic conductivity k F by K and named K the capillary conductivity depending on the moisture content, what might be expressed by the relative hydraulic conductivities through where s L = n L /n F is the liquid saturation. Furthermore, k L is the hydraulic conductivity of the pore liquid under fully saturated conditions, while κ L r addresses an empirical parameter that was given as in (22) 2 by Brooks and Corey [10] as a function of s L and a pore size distribution or wettability parameter λ.
Combining (21) with (20) yields the Richards equation 3 The classical equations in the light of the TPM

Theory of Porous Media in brief
As the classical equations described above, the equations governing porous media mechanics are based on a macroscopic, continuum mechanical description, where the porous solid material and the pore content are smeared out over the control volume. This leads to a continuum-mechanical setting, where the balance equations of the overall aggregate as well as of the individual components are given in the same domain, the latter coupled to each other by so-called interaction or production terms, respectively, compare Bowen [11,12] and Ehlers [13,14]. The mathematical basis of these equations is the Theory of Mixtures, compare Bowen [15], that has been extended by the concept of volume fractions. This concept goes back to the very early days of the investigation of porous and multi-component media and has been found and described by Woltman [16] and Delesse [17], also compare the short historical review of the TPM by the author [18]. In the following sections, a fluid-saturated porous material is discussed and compared to Darcy's law and the extensions of Forchheimer and Brinkman in the framework of a biphasic medium of solid and liquid, while the Richards equation describing a partially saturated material needs the investigation of a triphasic medium of solid, liquid and gas.

Biphasic materials
The biphasic medium ϕ under discussion is understood as a deformable porous solid material ϕ S saturated by a single pore fluid ϕ F . Both constituents, solid and fluid, are taken as materially incompressible components, meaning that their intrinsic densities ρ α R remain constant under isothermal conditions. For the solid skeleton, this means that although the intrinsic solid density remains constant, the partial solid density ρ S = n S ρ S R can vary through variations in the solid volume fraction. This model is composed by From this composition, volume fractions n α can formally be defined by the fraction of a volume element dv α of ϕ α over the total volume element dv of ϕ. Thus, This definition naturally includes the sum of all volume fractions yielding the so-called saturation condition α n α = 1.
As the following biphasic setting is mainly in line with Sect. 3.1 of Ehlers and Wagner [19], the governing set of balance equations is given as there by the mass and momentum balances of both constituents, solid and pore fluid via In these equations, ρ α is the partial density of ϕ α given as the product of the intrinsic (respectively, effective or real) density ρ α R and the volume fraction n α , T α is the partial Cauchy stress, g the gravitation vector and p α the so-called direct momentum production [13] coupling the momentum balances of solid and pore fluid.
Since the TPM provides a continuum-mechanical view onto porous-media problems, all terms have to be understood as local means representing the local averages of their microscopic counterparts. For example, the direct momentum production exhibits the volumetric average of the local contact forces acting at the local interfaces (pore walls) between the solid and the pore-fluid. In addition to the above, x α and x α are the local velocity and acceleration terms, while (·) α is the material time derivative following the motion of ϕ α . Assuming incompressible material properties in an isothermal setting by ρ α R = const., the mass balance (27) 1 reduces to the volume balance Combining the time derivative of the saturation condition n S + n F = 1 with the volume balances of solid and pore fluid leads to the saturation constraint of the overall aggregate or the incompressibility condition of materially incompressible solid and fluid materials reading where w F := x F − x S is defined as the seepage velocity. Note that the seepage velocity w F in deformable porous media and the filter velocity v F V are related to each other by the fluid volume fraction n F . Thus, Furthermore, note again that the fluid volume fraction n F is equivalent to the porosity, as the pore fluid in a fluid-saturated medium covers the whole pore space. As n F = 1 − n S and n S = n S 0 (det F S ) −1 with n S 0 as the solid volume fraction in the reference configuration at t 0 and det F S as the determinant of the solid deformation gradient F S , n F is determined through the solid motion via In continuum mechanics, the balance equations alone are not sufficient for the solution of initial boundary value problems (IBVP). Instead, the set of balance relations has to be closed by constitutive equations that have to fulfil the entropy inequality, which for isothermal processes reduces to the Clausius-Planck inequality where L α is the velocity gradient of ϕ α and ψ α the Helmholtz free energy of the constituents per unit constituent mass.
Based on the principle of phase separation [20] stating that the free energies of immiscible components of the overall aggregate, as on the microscale, do only depend on their own constitutive variables, the arrangement of an elastic solid and an incompressible fluid in an isothermal setting implies ψ S = ψ S (F S ) and ψ F = ψ F (−).
Multiplying the incompressibility constraint (29) with a Lagrange multiplier Λ ≡ p F R and adding the result to (32) yield the entropy inequality in its final shape: Therein, L α has been substituted by its symmetric part D α on the basis of the symmetry assumption for T α .
are the so-called extra terms of stress and direct momentum production withp F +p S = 0. The quantity I is the second-order identity tensor, while p F R is the excess pore pressure meaning the pressure exceeding the atmospheric pressure that, in short, is simply called the pore pressure. Note in passing that the solid extra stress coincides with the so-called effective stress T S eff defined as that part of the overall stress state that produces solid deformations by convenient constitutive equations, compare Ehlers [21].
Exploiting the entropy inequality yields by the use of the Coleman-Noll procedure [22] that the following results hold: Therein, the equation for T S E = T S eff has been obtained from the equilibrium part of (33), where D S can be arbitrary and D F and w F vanish, while the extra or frictional terms T F E = T F fric andp F E =p F fric are found from the dissipational part as functions of D F and w F , respectively, with positive definite forth-and second-order constitutive tensors 4 D v and S v given as with ( I ⊗ I ) 23 T as the fourth-order fundamental or identity tensor. Furthermore, note again that μ F = n F μ F R is the partial shear viscosity of the pore fluid, while different possibilities to describe permeability properties alternatively to (35) 2 can be taken from (5). Finally, note in passing that a scalar-valued hydraulic conductivity like k F is only valid in case of an isotropic pore structure, while anisotropic pore structures require a symmetric tensor K F of hydraulic conductivities that reduces in case of isotropy to k F I.
Concerning the solid skeleton, a geometrically linear version of the elasticity law satisfying (34) 1 reads where the extra Cauchy stress T S E or the effective stress T S eff , respectively, approximately coincides with the linearised extra or effective stress σ S eff . Furthermore, ε S = 1 2 (F S + F T S ) − I is the linearised solid strain tensor, and μ S and λ S are the Lamé constants of the porous solid.
Given the above constitutive equations for the solid and fluid extra stresses and the extra momentum supply, the governing equations of the binary model result in Therein, (37) 1 is a modification of the overall volume balance (29), while (37) 2 and (37) 3 are results of (27) This yields the so-called modified Eulerian description, where the pore fluid is not described with respect to a fixed frame but with respect to the frame of the moving skeleton. Applying the modified Eulerian description also to the fluid volume balance, (28) additionally yields Furthermore, as the external loadt acting at the Neumann boundary of the saturated porous solid material cannot be split into solid and fluid portions, the solid and fluid momentum balances have to be added. Thus, (37) 2 is substituted with the aid of (38) by with ρ := α n α ρ α R as the so-called mixture density of the overall aggregate. Note in passing that under quasi-static and creeping-flow conditions with (u S ) S ≈ 0 and (v F ) F ≈ 0, the left-hand side of (40) vanishes and the overall momentum balance reads These equations are sufficient to conclude to the classical fluid flow descriptions of Darcy, Brinkman and Forchheimer.

The Darcy and Darcy-Brinkman equations
Concerning the reduction in the setting (37) towards the findings by Darcy and the modification by Brinkman, one has to be aware that the classical equations proceed from a non-deforming solid skeleton, such that u S ≡ 0. Thus, (u S ) S is dropped from (37) 1 , while (37) 2 does not play any role and is dropped completely. Thus, only remains. Furthermore, as the solid is non-deforming or rigid, n S is constant in time and additionally constant in x, once the solid material has a homogeneous pore structure. Assuming these conditions, (42) can be divided by n F yielding where w F = v F together with x S = (u S ) S ≡ 0 has been used. Based on (35), the frictional fluid stress T F fric and the frictional momentum productionp F fric or the drag force read Combining (43) and (44) results in where div grad T v F = grad div v F = 0 has been used on the basis of (43) 1 . In (43) 2 the frictional force div Following this, it is worth to check this competition by the use of a dimensional analysis, cf. Ehlers et al. [23].
Introducing dimensionless quantities by * (·) such that * with L as the characteristic length and V as the characteristic velocity yields the following relations between dimensionless and non-dimensionless gradient and divergence operators: Following this, the frictional force and the drag force can be reformulated as Given these equations, the quotient of the prefactors of the frictional and drag forces is a measure for the importance of these terms in different versions of the fluid momentum balance. Thus, one concludes to where (5) has been used. Once K S is known, the relation between the frictional force and the drag force is mainly governed by the characteristic length L, as the viscosity μ F R drops out. Understanding K S as the internal length parameter of the porous microstructure, investigations at the microstructural domain with L √ K S yield the frictional force to be dominant, while investigations at the macroscopic domain with L √ K S yield the result that the frictional force is negligible compared to the drag force. However, the frictional force can always be taken as a corrector to the drag force, compare Auriault [5].
To set an example, consider a geomaterial like a dense sand with a porosity of n F = 0.35 and an intrinsic permeability of K S = 2 × 10 −11 m 2 , while L is set to 0,1 m, meaning, for example, a typical width or height of an experimental setup. In this case, K S over n F L 2 yields 5.71 × 10 −9 .
Thus, the frictional force is negligible and (45) reduces to Darcy's filter law In the subsurface, one usually finds creeping flow conditions, such that the acceleration (v F ) F is negligible in comparison with the other terms in (50). Following this, (50) can be solved with respect to the filter velocity v F V = n F v F yielding Darcy's law, compare (4): In conclusion, Darcy's filter law is recovered from the TPM by the application of a binary model of one porous solid and one fluid. Furthermore, both constituents have been assumed to be materially incompressible with constant ρ α R in space and time. This led to Eq. (37) that could further be reduced by the assumption of a non-deforming or a rigid skeleton, respectively, where n F is constant in space and time, compare (43). In addition to this, it could be shown by a dimensional analysis that the frictional force div T F fric can be neglected compared to the drag forcep F fric . Finally, the inertia term in the fluid momentum balance (50) dropped out as a result of the creeping-flow assumption, such that (50) could be solved with respect to the filter velocity v F V .
Remark As Darcy's law (51) has been obtained without the need to use the assumption of a non-deforming skeleton, Darcy's law also holds without further restrictions when v F V is identified as n F w F including the skeleton motion: The Darcy-Brinkman equation In contrast to the Darcy law valid at the macroscopic range with L √ K , investigations at the microstructural domain proceed from L √ K , such that the drag force is negligible compared to the frictional force. In this case, the fluid momentum balance reduces to the Navier-Stokes equation, where also the acceleration terms might have to be considered. However, in between these borderline cases, also L ≈ √ K is possible, such that both the frictional and the drag force have to be taken into consideration. In this case and under creeping flow conditions with (v F ) F ≈ 0, (45) reads Solving this equation with respect to the pressure gradient gives Assuming that the body forces do not matter, this equation yields with the aid of (5) that thus recovering the Brinkman equation (19).

Remark 1
In contrast to the Darcy equation, the Brinkman equation cannot be adopted to fluid flow in deformable porous solids without discussion. This is due to the fact that the term div T F fric = μ F R Δ v F was found on the basis of (43) 1 . However, for deformable solid skeletons, (43) 1 has to be substituted by (42) 1 , meaning that div v F does not necessarily vanish, such that div In other words, using the Brinkman equation for flow in deformable porous solids includes the additional assumption that μ F R grad div v F approximately vanishes.
Remark 2 As has been stated after (19), Brinkman wanted to show that his equation containing a corrector term compared to the original Darcy equation was appropriate to describe situations with low and high permeabilities. This usually describes a situation, where a transition from a porous domain towards a free flow domain is considered. On the other hand, the above dimensional analysis also underlines that, depending on the size of the pore structure, the Navier-Stokes term in (55) becomes apparent when microstructural domains are of interest.

The Darcy-Forchheimer equation
For the validity of the Forchheimer equation, the same conditions are valid as for the Darcy equation. In particular, the viscous force is also negligible compared to the drag force. Neglecting gravitational forces, the Darcy law (51) can be written as where (5) has been used. Now, (56) is an equation for the pressure loss as a linear function of the filter velocity. However, in cases where this linearity is not found in the experimental evidence, a nonlinear pressure loss can be modelled by a variation in (35) 2 . Maintaining the thermodynamic validity of the model, (34) 3 can be extended towardsp In addition to (35) 2 , S v and S q can be found as positive definite tensors where the nonlinear velocity term is controlled by the tortuosity parameter B S describing the inherent irregularity of the pore channels. Neglecting inertia and gravitational terms and disregarding the viscous force compared to the drag force, an insertion of (57) and (58) in (43) 2 yields the Darcy-Forchheimer equation (9) that has been expressed in terms of the filter velocity v F V = n F w F and grad p F R instead of grad h: Proceeding as in Sect. 2.2, |v F V | has to be expressed by |grad p F R |, such that (59) finally reads This equation is identical to the Forchheimer equation (15). Finally, note in passing that Hassanizadeh and Gray [24] came to a comparable result by the use of volume averaging techniques within the hybrid mixture theory.
Remark In contrast to the Darcy equation, the Darcy-Forchheimer equation proceeds, on the one hand, from negligible body forces and, on the other hand, from the addition of a quadratic filter-velocity term. As the viscous force has not been included, the Darcy-Forchheimer equation can also be applied to deformable media, whenever the limiting restrictions necessary to obtain Darcy's law are valid.

Triphasic materials
For the investigation of the continuum-mechanical validity of the Richards equation, the consideration of a triphasic medium of solid skeleton, pore water and pore gas is necessary. According to Sect. 3.2, the triphasic medium ϕ under discussion is understood as a deformable porous solid material ϕ S that is saturated by two immiscible pore fluids, a materially incompressible pore liquid ϕ L and a compressible pore gas ϕ G . Thus, the model is composed of The balance equations for mass and linear momentum can directly be taken from (27) reading where the only difference between (27) and (62) where, in addition to the volume fractions n β , s β characterises the saturations of ϕ β as fractions of the porosity n F . In a general description including a deformable solid, such that n F = n F 0 with n F 0 as the reference porosity at t 0 , only n F can be computed as a function of the solid deformation gradient F S through n F = 1 − n S with n S = n S 0 (det F S ) −1 , while the liquid and gas saturations have to be found by the use of (62) 1 through As the pore liquid is assumed to behave materially incompressible at constant temperature, the first equation of (64), after having been divided by ρ L R , reduces to a volume balance, while the second equation reflects that the rate of the partial density ρ G of a compressible gas includes the rates of both the volume fraction n G and the intrinsic density ρ G R . At this point, one has to discuss basic features of unsaturated soil models. Several authors proceed from a so-called static gas phase and argue that this is not only easier to handle compared to the fully triphasic description but also sufficiently exact, compare, for example, Callari and Abati [25] or Moldovan et al. [26]. However, using the static gas-phase model means that the gas pressure p G R is always set to the atmospheric pressure. As p G R and ρ G R are usually coupled, for example, by the ideal gas law, this assumption furthermore includes ρ G R = const. at constant temperature, such that also (64) 2 reduces to a volume balance. But the assumption of a static gas phase includes further features. Taking p G R as the excess gas pressure, such that p G R ≡ 0, the capillary pressure p C = p G R − p L R reduces to p C = −p L R . In standard drainage and imbibition processes, p C is defined as the positive difference between the non-wetting and the wetting fluid, gas and liquid. Thus, p L R must be negative describing the capillary suction. Furthermore, as equation (64)   is used in numerical models to find the values of p G R , this equation is no longer needed and the triphasic description reduces to a quasi-triphasic or a biphasic one, respectively. As a result of the above, the effect of capillary suction, where the excess pressure of the pore gas is negative, cannot be modelled.
To set an example of this effect, consider the Liakopoulos test (Liakopoulos [27]) that has been investigated by the ALERT group 1 as a benchmark (e. g. Klubertanz et al. [28]) and was modelled and computed, for example, by Ehlers et al. [29], compare Figs. 3, 4 and 5 taken from that article. In these figures, it is clearly seen for the example of the leakage of pore fluid from a fully saturated elastically deformable soil column of very fine Del Monte sand under the effect of gravitational forces that the difference between the assumption of a static gas phase and a fully triphasic medium matters a lot. From Fig. 3, one observes that the capillarity effect combined with the drag force retains the pore water much longer as the drag force alone. Figure 4 furthermore exhibits the capillarity effect through the excess gas pressure p G R that goes down to around − 6,200 Pa at approximately 30 min after having opened the outflow, while, in the static gas case, p G R remains zero.
Finally, the advantage of the fully triphasic model compared to the biphasic simplification is shown in Fig. 5, where the water efflux is plotted versus time. In particular, a comparison of the experimental data by Liakopoulos with the bi-and triphasic formulations reveals that the biphasic computation with the static gas phase overestimates the water efflux at the beginning of the draining process, whereas the fully triphasic computation yields a slight underestimation of the primary outflow. In contrast to the experimental evidence, the biphasic formulation predicts the final distribution of the liquid saturation at t = 5 h, while the triphasic formulation fits the experiments very well. In this formulation, the final saturation distribution is reached after approximately 88 h, compare Fig. 3.
From these results, it can easily be concluded that things are as complicated as they are, but not easier. Making them easier always yields a loss of information that might lead to wrong results and, in conclusion, to wrong engineering decisions.
As the triphasic setting is much more complicated than a biphasic one but mainly in line with Sect. 5.2 of Ehlers [14] and Sects. 2.3-2.5 of Ehlers [21], these sections can be considered for further details. In the triphasic description, the saturation condition is obtained in analogy to (29) through the combination of the time derivative of the saturation constraint n S + n L + n G = 1 with the volume and mass balances (64) leading to 0 = n S div x S + n L div x L + grad n L · w L +  The constitutive equations of the triphasic model have to fulfil the same Clausius-Planck inequality (32) as those of the biphasic model with the only difference that three components instead of two have to be considered: Based on the principle of phase separation [20], the free energies ψ α are concluded to only depend on their own constitutive variables. As a result, ψ S depends as in the biphasic model on F S , while ψ G is a function of the intrinsic gas density ρ G R . In analogy to the biphasic model, ψ L does not depend on ρ L R as ψ F did not depend on ρ F R . However, ψ L depends on s L . Thus, such that While ψ L depends on s L , s L depends on n L and n F = 1 − n S 0 (det F S ) −1 . Therewith, the liquid saturation is coupled to the solid deformation what has to be taken into consideration by where (69) has been obtained by the use of the mass and volume balances.
Multiplying the saturation condition (65) with a Lagrange multiplier Λ and adding the result to (66) yield the final form of the Clausius-Planck inequality reading where (68) and (69) have been used together with Exploiting the entropy inequality can easily be done in two steps following the Coleman-Noll strategy (Coleman and Noll [22]). In the first step, the thermodynamic equilibrium is considered, meaning that (70) has to vanish for arbitrary values of D S , D L and D G , and additionally for arbitrary values for (ρ G R ) G , w L and w G . As a result, one concludes to the equilibrium terms As (72) 3 is correct, whenever the term in parentheses equals p L R , one concludes from with p C as the capillary pressure. Given (72) and (73), one furthermore finds Taking a closer look at the term in parentheses of (74) 2 , it is easily concluded that recovers Dalton's law [30]. With (74) 2 in mind, it is concluded that p F R = p S R , as the pore pressure p F R acts on both the overall pore content and the porous solid. With Dalton's law, (74) 2 reads as the effective stress. As in the biphasic model, T S eff is compatible with the Hookean law in the framework of linear elasticity, compare (34) 1 and (36). Finally,p Once the equilibrium terms are known, the non-equilibrium parts of T L , T G ,p L andp G have to be found from the dissipation inequality Following this, one concludes to Proceeding from (71), (72) 1 and (73) 2 , the proportionalities in (79) can be expressed as follows: where T β fric andp β fric are the frictional liquid and gas stresses and drag forces that are defined as in (35) through 4 D β v = 2μ β ( I ⊗ I ) 23 T , However, in contrast to (35) 2 , (81) 2 is not governed by hydraulic conductivities k β but by relative hydraulic conductivities k β r = κ β r k β also known as capillary conductivities. These terms are related to the hydraulic conductivities k β of ϕ β through the dimensionless relative permeability factors κ β r that depend on the degree of saturation, cf. (22). Finally, note in passing that a Forchheimer addition of a quadratic term as in (57) and (58) is possible but has not been taken into consideration here.
In the general non-equilibrium case, the above equations finally result in Inserting these equations into the momentum balances (62) 2 for different components of the overall model, one obtains for the solid skeleton wherep S = −p L −p G together with (82) 3,4 has been used. Analogously, the pore-fluid equations can be found by the use of Dalton's law (75) as Note again that in a numerical scheme proceeding from a modified Eulerian description, the liquid and gas accelerations included in (84) have to be modified via also compare (38).
Applying the above full triphasic model to numerical computations by use of the Finite Element Analysis (FEA), the primary variables are the solid deformation u S obtained from the weak form of the solid momentum balance (83), the fluid velocities v β obtained from the weak forms of the fluid momentum balances (84) and the effective fluid pressures p β R obtained from the weak forms of the liquid volume balance and the gas mass balance given in (64). As for the biphasic model, the solid momentum balance is usually substituted by the overall momentum balance of solid, liquid and gas. This is again due to the fact that it is generally not possible to split the external loadt acting at the Neumann boundary into portions acting separately on the solid, the liquid and the gas. Thus, by addition of (83) and (84), the overall momentum balance reads where (85) has been used. Beneath the constitutive equations for the effective stresses T S eff as the mechanical part of the solid stresses T S , the fluid frictional stresses T β fric and momentum productionsp β fric are given to close the system of governing equations. However, as only n S can be determined by F S = I + Grad S u S with Grad S (·) = ∂(·)/∂x 0 and n F = 1 − n S , the fluid volume fractions n β = s β n F or the saturation s β , respectively, must be found from an additional constitutive equation, for example, from the Brooks and Corey law [10] as a function of the capillary pressure p C .

Simplifications of the full triphasic model
The above dynamic triphasic model can be used for various numerical applications and can easily be extended, for example, towards finite solid deformations and towards nonlinear elastic, viscous and plastic solid material properties. These applications usually proceed from a quasi-static setting, where the acceleration terms in the above momentum balances are neglected. On the other hand, when earthquake problems come into play, the solid skeleton vibrates so quickly that the fluid has no possibility to creep away, such that the fluid accelerations are set approximately identical to the solid vibrations.
However, three different basic loading scenarios of porous media can be distinguished, those, where the loaded solid skeleton drives the fluid flow as in consolidation or base failure problems, or those, where the fluid flow drives the solid deformation as in hydraulic fracturing. The third group is none of the above, but it is in the centre of the following simplifications. Here, the solid skeleton is not loaded except of gravitational forces and fluid flow resultants, more precisely fluid pressures, and these effects are not large enough to induce solid deformations. For the present article, this last group of problems is essential for the investigation of the Richards equation. Following this, the full triphasic model will be reduced in 4 steps.
Step (1): Quasi-static scenarios In these scenarios, inertia forces can be neglected compared to the other forces included in the momentum balances. This means the solid skeleton is assumed to be only under static loading and the fluids exhibit creeping flow conditions. Thus, the momentum balances (84) and (86) of the full model reduce to Remark This model built from the above momentum balances and the volume and mass balances (64) still needs the full set of primary variables as the frictional fluid stresses still proceed from D β through the gradient of v β .
Step (2): Neglect of frictional fluid forces Once the frictional forces div T β fric can be neglected compared to the drag forcesp β fric , (87) simplifies to As there are no frictional forces, the primary variables v β are no longer needed and (88) 2 and (88) 3 can be solved with respect to the filter velocities n β w β . Thus, where (80) 2,3 and (81) 2 have been used. This formulation, given by (88) 1 and (89), is known as the u-pp formulation with the primary variable u S , p L R and p G R . In a numerical setting, the overall momentum balance (88) 1 can be used for the determination of the solid displacement u S , while (89) 1 and (89) 2 are inserted in the liquid volume balance and the gas mass balance, compare (64), yielding with filter velocities n β w β from (89).
Step (3): Static gas phase Proceeding from a static gas phase induces the assumption of a constant gas pressure equivalent to the ambient or the atmospheric pressure, respectively. As a consequence, the intrinsic gas density ρ G R is constant under isothermal conditions and the gas mass balance reduces to a volume balance reading In a numerical setting, this equation could be used to compute the gas pressure p G R . However, as the gas pressure is constant, this equation is not needed anymore and the set of governing equations reduces to a quasi-triphasic setting yielding As all pressures have been normalised at the atmospheric pressure p 0 , such that p β R = 0 at p 0 , p C = −p L R .
In the framework of the usual p C over s L curves describing standard drainage and imbibition processes like the Brooks and Corey [10] or the van Genuchten [31] laws, p C must always be positive or zero, such that p L R can only describe a suction with p L R ≤ 0, while, on the other hand, spontaneous drainage and forced imbibition curves also exhibit negative values of p C , compare Blunt [32]. Finally, note in passing that the term in brackets governed by p L R /n L has been found from the Clausius-Planck inequality as a thermodynamical restriction of the triphasic model. However, as many authors working on triphasic media do not check thermodynamical restrictions, this term is usually ignored.
Remark Equations (92), although obtained for a triphasic medium with static gas phase, provide a hindsight towards biphasic media. Setting s L = 1 implies n L ≡ n F , such that the first two equations of (92) correspond to (39) and (40). Furthermore, as s L = 1 also includes s G = 0, the last term of (92) 3 governed by p L R /n L vanishes, such that (92) 3 reduces to Darcy's law with k L r ≡ k F , compare (22) and (52).
Step (4): Rigid solid skeleton The rigid-skeleton assumption is based on a non-deforming solid matrix, such that u S = 0. As a result F S = I and det F S = 1. This yields a constant porosity n F = n F 0 = const., such that n L = s L n F is only governed by the liquid saturation. As a further consequence of the rigid skeleton assumption, equation (92) 1 governing u S can be dropped and (92) reduces to Note that while stepping from (92) 3 to (93) 2 , the term governed by p L R /n L has been neglected in view of the Richards equation, such that (93) 2 equals Darcy's filter law, however, applied here to the liquid phase of a triphasic material, compare (52). It should furthermore be noted that, as a further result of u S = 0, (n L ) S + n S div (u S ) S reduces to ∂n L /∂t.

The Richards equation
Based on (93), the Richards equation can easily be found. The assumptions of creeping flow conditions in a partially saturated rigid solid skeleton with a static gas phase are applied. Furthermore, as γ L R = const., (93) 2 can be written as compare equation (2). Inserting (94) in (93) 1 yields thus recovering the Richards equation. Note in passing that even in case of a homogeneous pore-size distribution with constant k L , k L r is generally not constant as it depends on x through s L , compare (22).

Conclusions
In hydraulic engineering as well as in geomechanics, the well-known equations of Darcy, Forchheimer, Brinkman and Richards are widely used as constitutive equations without considering their validity range. Following this, it was the goal of the present article to show that these equations perfectly fit into the framework of the Theory of Porous Media, although they have been found only by empirical analyses and mathematical analogies. On the other hand, the TPM as a general tool for the description of multi-component and multiphasic materials offers a wide basis for the solution of boundary-value problems, especially in the fields of geomechanics and biomechanics and in all other fields, where coupled problems of porous solids with fluid pore content have to be considered. This naturally includes the possibility to distinguish between miscible and immiscible pore fluids.
In the field of fully water-saturated solids, for example, soils, the TPM provides a quite simple basis for the description of coupled solid-fluid problems and furthermore guarantees the thermodynamic validity of the elaborated models, no matter whether or not the solid is deformable or if the fluid (liquid or gas) is an inert fluid or a fluid mixture. As has been shown, the Darcy, Forchheimer and Brinkman equations perfectly fit into the character of the biphasic TPM. As Darcy, Forchheimer and Brinkman had only non-deforming solid skeletons in mind when they presented their famous equations, the derivation of these equations followed their assumption. Further restrictions are the assumption of homogeneous and isotropic pore structures, such that the hydraulic conductivity could be assumed a single constant scalar. As the Darcy and Forchheimer equations are usually applied to macroscopic problems, it could furthermore be shown by a dimensional analysis that frictional forces could be neglected in comparison with drag forces. While the Darcy equation gives a linear relation between the filter velocity and the gradient of the hydraulic head, the Forchheimer equation couples the gradient of the hydraulic head with a quadratic function of the filter velocity, thus including nonlinearities into the description. In contrast to Darcy and Forchheimer, Brinkman did not only include the drag force but also the frictional force in his equation and therewith presents both a possibility to include a transition from porous media flow towards free flow conditions, on the one hand, and to describe microstructural effects in the pore fluid flow relations, on the other hand. To set an example, consider hydraulic fracturing scenarios, where a fracking fluid is pressed into a porous soil or rock with the goal to drive a fracture. In this case, a Darcy-type pore fluid flow appears in the porous matrix with transition towards a Navier-Stokes-type flow in the fracture, compare, Ehlers and Luo [33,34]. Finally, it could be shown that the Darcy and Forchheimer equations can also be applied to deforming solid skeletons in the realm of their basic restrictions, while the Brinkman equation needs an additional assumption.
While the Darcy, Forchheimer and Brinkman equations could easily be found from biphasic TPM models, the Richards equation describes a partially saturated soil, where the pores are filled with two immiscible fluids, the pore liquid and the pore gas.
Proving the validity range of the Richards equation made it necessary to firstly investigate a general triphasic model of solid, liquid and gas that had to be broken down in 4 steps to obtain Richards' equation. These steps consist of (1) quasi-static scenarios, (2) the neglect of frictional forces, (3) a static gas phase and (4) a rigid solid skeleton.
The above statements make clear that taking the classical hydromechanical equations as constitutive laws does not necessarily lead to satisfactory solutions of related boundary-value problems, unless the range of validity of the classical laws is not fully considered, also compare the work by Gray and Miller [35].
Funding As Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.