Generalized Buckley–Leverett theory for two-phase flow in porous media

Hysteresis and fluid entrapment pose unresolved problems for the theory of flow in porous media. A generalized macroscopic mixture theory for immiscible two-phase displacement in porous media (Hilfer 2006b Phys. Rev. E 73 016307) has introduced percolating and nonpercolating phases. It is studied here in an analytically tractable hyperbolic limit. In this limit a fractional flow formulation exists, that resembles the traditional theory. The Riemann problem is solved analytically in one dimension by the method of characteristics. Initial and boundary value problems exhibit shocks and rarefaction waves similar to the traditional Buckley–Leverett theory. However, contrary to the traditional theory, the generalized theory permits simultaneous drainage and imbibition processes. Displacement processes involving flow reversal are equally allowed. Shock fronts and rarefaction waves in both directions in the percolating and the nonpercolating fluids are found, which can be compared directly to experiment.


Introduction
A predictive theory of two-phase flow in porous media is a longstanding challenge for theoretical physics (Muskat 1937, Kozeny 1953, Scheidegger 1957, de Wiest 1969, Bear 1972, Kadanoff 1985, de Gennes 1988, Lenormand et al 1988. Displacement processes of two immiscible fluids inside a rigid porous medium have motivated and advanced several subfields of theoretical physics ranging from percolation theory (Kirkpatrick 1973, Chandler et al 1982, Stauffer and Aharony 1992, Sahimi 1993) and pattern formation (Weitz et al 1987, Oxaal 1991, Ferer et al 2003, Babchin et al 2008 to wetting phenomena (de Gennes 1985, Bonn et al 2009 or nonlinear field theories (Hopf 1969,Glimm et al 1982, LeFloch 2002. Macroscopic models for two-phase flow in porous media on the scale of 10 2 -10 3 m (e.g. petroleum reservoirs) generally start from the well known Buckley-Leverett (BL) equation (Buckley and Leverett 1942). Although this nonlinear partial differential equation is nowadays universally accepted as a starting point in theoretical physics (Bogdanov et al 2003), mathematics (Bourgeat 1997), numerics (Chen et al 2006) and engineering (Lake 1989), it has never been derived rigorously from the underlying microscale equations (Hilfer 1996). Its solutions cannot reproduce simple experimental observations such as capillary desaturation or hysteresis, 3 because it is valid only in the limit of vanishing flow velocities when fluid-fluid interfaces or contact lines do not move (Hilfer 2006a). Our objective in this paper is to deduce and solve generalized macroscopic equations valid for finite velocities, which include the traditional BL theory in the limit of vanishing velocities. Realistic flows with finite velocities produce motion of fluid-fluid interfaces and contact lines (Lenormand et al 1988). Experimentally this gives rise to macroscopic hysteresis between imbibition and drainage processes and to spatial and temporal variations of residual, trapped and immobile fluid regions (Payatakes 1982). Mobilization of these trapped residual fluids is also technologically very important.
Distinguishing and identifying the physical importance of percolating (= hydraulically connected) and nonpercolating (= hydraulically isolated) fluid parts (Hilfer 1998, Hilfer andBesserer 2000) within the classical BL problem is the main objective, physical concept and innovation of this paper. Explicit generalized equations will be extracted from the equations introduced in Hilfer (2006aHilfer ( , 2006bHilfer ( , 2006c) using a hyperbolic approximation. It will then be possible to formulate, study and solve the analogue of the classical BL frontal advance problem, and to introduce and solve in addition a frontal advance problem with subsequent flow reversal.
Given these objectives we recall that the classical BL problem provides a simple example where the nonlinear partial differential equations can be solved analytically without numerical approximations (Buckley and Leverett 1942, Aziz and Settari 1979, LeVeque 1992. Let us also recall that the numerical solutions are difficult to obtain, because classical solutions of the BL problem may develop discontinuities from smooth initial data, as it is characteristic for nonlinear hyperbolic conservation laws (Lax 1957, Glimm et al 1982. One needs and utilizes the analytical BL solution therefore routinely to test numerical algorithms (Aziz and Settari 1979, LeVeque 1992, Chen et al 2006. Reformulating the equations in Hilfer (2006aHilfer ( , 2006bHilfer ( , 2006c with the equal pressure assumption and other experimentally justifiable assumptions yields a fractional flow formulation analogous to the traditional BL theory. It is found that shock waves may appear from this formulation as in the BL case. A difficulty, which is not present in the traditional BL theory, arises from spatially discontinuous flux functions due to hysteresis. Monotone flux functions will be discussed in this paper, nonmonotone functions in a forthcoming publication. The theory of immiscible two-phase flow in porous media is a challenging and longstanding open problem of theoretical physics, because of the multiple scales and the nonlinearity of the governing equations. On microscopic scales wetting properties and roughness of surfaces are important. On mesoscopic scales the flow is determined by the complex pore structure and mobile fluid-fluid interfaces. On larger scales there is considerable uncertainty about the important physical observables, as well as geometrical and transport parameters that govern the properties of the permeable fluid-solid mixture. As a result the linear response (or mean field) theory for single-phase flow in porous media is well established to be Darcy's law, while the linear response theory for two-phase flow is still unknown at present. The present paper introduces important physics insights and linear response concepts to overcome this impediment. This paper is structured as follows. In section 2 a brief review and the fractional flow formulation are given. In section 3 pertinent approximations are formulated. In section 4 the number of equations is reduced. The general solution is obtained based on the method of characteristics in section 5. The Riemann problem is formulated in section 6. It is solved and discussed in section 7, and generalizes the BL problem (Buckley and Leverett 1942) in the sense of accounting for percolating and nonpercolating fluid phases. The paper concludes with several 4 remarks in section 8. For the convenience of the general reader, the widely used traditional theory (BL theory) is recapitulated in an appendix, which highlights the analogy between the traditional approach and this paper.

Balance laws and constitutive assumptions
In this section a brief description of the model, which accounts for percolating and nonpercolating fluid phases, is given. For details we refer to preceding works. The mathematical model was originally formulated in Hilfer (2006aHilfer ( , 2006bHilfer ( , 2006c. It is based on ideas introduced earlier in Hilfer (1998) and Hilfer and Besserer (2000). A self consistent closure condition was proposed and used in numerical simulations in Hilfer and Doster (2010) and Doster et al (2010). In this paper the model is restricted to one spatial dimension.
Percolating and nonpercolating fluid parts are distinguished as separate phases and, hence, a two-phase system is treated as a four phase system (see Hilfer (1998), Hilfer and Besserer (2000), Hilfer (2006aHilfer ( , 2006b for definitions and details). Following the notation of Hilfer (2006b) we shall name one fluid water, and identify it with index W, and the other fluid oil, and identify it with index O. The percolating water is identified with index i = 1, the nonpercolating water with i = 2, the percolating oil with i = 3 and the nonpercolating oil with i = 4. For an incompressible porous medium volume conservation requires S 1 + S 2 + S 3 + S 4 = 1 (1) to hold, where S i (x, t) is the saturation of phase i as function of position x and time t. The mass balance for fluid phase i ∈ {1, 2, 3, 4} in a porous medium is expressed in differential form as where i (x, t), φ(x, t), v i (x, t) are mass density, porosity and velocity of phase i as functions of position x ∈ [x , x r ] and time t ∈ R + and M i is a source or sink term for the corresponding phase. For incompressible fluids and rigid homogeneous media the porosity and the densities do not depend on position x and time t for i = 1, 2, 3, 4. One has where W , O are the densities of water (wetting fluid) and oil (non-wetting fluid). The source terms M i represent mass exchange between percolating and nonpercolating phases through breakup and coalescence. They are assumed as where η 2 , η 4 are constants and the shorthand ∂ t = ∂/∂t is used. The parameters S W * , S * 2 , S * 4 are defined by where S W dr , S O im are limiting saturations for S 2 , S 4 and (x) denotes the Heaviside unit step function. The momentum balance is written for phase i ∈ {1, 2, 3, 4} as where D i /Dt = ∂/∂t + v i ∂/∂ x denotes the material derivative for phase i, i the stress tensor in the ith phase, F i the body force per unit volume acting on the ith phase and m i denotes the momentum transfer into phase i from all the other phases. The material derivatives are neglected here because realistic subsurface flows have very small Reynolds numbers (Hilfer 1996). The stress tensors for the four phases are specified as where P 1 (x, t) and P 3 (x, t) are the fluid pressures in the percolating phases. The constants P * 2 , P * 4 and exponents γ , δ are constitutive parameters, which account for the energy density stored in the common interface with the surrounding percolating phase. The body forces are given by gravity and capillarity. They are specified as with constitutive parameters a , b , α, β > 0. The capillary terms in F 2 , F 4 are body forces derived from a capillary potential (see Hilfer 2006aHilfer , 2006bHilfer , 2006c. They serve to retain nonpercolating phases inside the medium. The angle 0 ϑ π/2 is the angle between the direction of the column and the direction of gravity with ϑ = π/2 corresponding to alignment. Finally, the momentum transfer terms are assumed to be given by linear viscous drag characterized by constitutive resistance coefficients R i j through the equations where R 12 = 0 and R 34 = 0 was used, because there is no common interface and hence no direct viscous interaction between these phase pairs. The index i = 5 represents the rock matrix. For physical motivation of the constitutive assumptions we refer the reader to the original publications (Hilfer 2006a(Hilfer , 2006b(Hilfer , 2006c.

Reformulation of the model
The balance laws (1), (2) and (7) are rearranged in terms of volume flux densities. The righthand side of the momentum balance (7) is linear in the velocities v i and can be written as  where the components of the generalized resistance matrixR arẽ and the shorthand notation was used. Inserting equations (12) into (7), solving the set of equations for the linearly appearing velocities, and multiplying with the corresponding phase volume density where phase volume flux densities Q i and the mobility matrix have been introduced. The components of are given as 7 in terms of the inverse ofR. The mass balances (2) are combined linearly and divided by the densities, and volume conservation (1) is applied to obtain where Q W = Q 1 + Q 2 and Q O = Q 3 + Q 4 . Note, that equation (17a) is identical to the transport equation of the traditional theory and (17d) to its elliptic pressure equation as formulated in many textbooks (Aziz and Settari 1979).

Hyperbolic approximation
In this section a hyperbolic limit of equation (17a) analogous to the BL theory is formulated.
The approximations lead to a fractional flow formulation in section 4.

Immobility of nonpercolating fluids
To ensure that the nonpercolating phases are immobile, it will now be assumed that for all pairs (i, j) with i ∈ {1, 2, 3, 4} and j ∈ {1, 2, 3, 4, 5} such that (i, j) = (2, 5) and (i, j) = (4, 5). It will also be assumed that holds. In the limit R 25 → ∞, R 45 → ∞ only the components 11 , 13 , 31 , 33 of the mobility matrix differ from zero. Hence fluxes are approximated as where the nonzero components of the mobility matrix are explicitly given as in terms of saturations. Here Onsager reciprocity R 31 = R 13 is assumed to hold for the viscous cross-coupling. The approximation of immobile nonpercolating phases is inspired by the residual decoupling approximation (Hilfer 2006c). It differs from the latter, because equations (20c) and (20d) follow from (15), and are now identically fulfilled. The approximation may be justified physically by the observation that the motion of contact lines on the internal surface requires capillary forces to be overcome, and this creates an additional resistance that is much higher than the viscous drag. Note that for R 13 = 0, a comparison of the mobilities in equation (21) with mobilities of the traditional theory yields Hilfer (2006bHilfer ( , 2006c) where k denotes the permeability and µ W , µ O the viscosities of water and oil respectively.

Viscous domination
It is assumed that viscous drag dominates the momentum transfer such that holds true. Then the resistance matrix becomes a constant matrix. The assumption is reasonable in many cases, because resistance coefficients are typically of order 10 8 kg m −3 s −1 while momentum transfer due to mass transfer is of order 10 3 kg m −3 s −1 . However, if S W → S W * or ∂ t S W → ∞, this assumption breaks down. Here we assume |S W * − S W | < 10 −4 and ∂ t S W < 1 s −1 . The assumption |S W * − S W | < 10 −4 is reasonable, because the saturations S W * and S W are rarely known with a precision better than 10 −4 . The assumption on ∂ t S W is justified, because the divergence resulting from shock fronts is smeared out in experiments.

Equal pressure assumption
In many applications like oil recovery the applied pressures are large enough to neglect differences in the pressures of the two percolating phases. Here we use this observation to assume 9 and hence 1 = 3 (26) for the stress tensors of the percolating phases. This is a strong assumption, which eliminates diffusive capillary terms. The same assumption is made in the traditional BL theory (see appendix equation (A.15)). The system becomes nonlinearly hyperbolic and shocks may occur even for smooth initial data.

Model reduction and fractional mobility functions
In this section a fractional flow formulation is derived. Fractional flow and mobility functions are identified in analogy to the BL theory (see appendix).
The approximations (18), (19), (23) and (25) are inserted into equations (17).  (20) and (25) into (27) where are the mobilities of water W and oil O. The fractional flow function f W and the fractional mobility λ are introduced in analogy to the traditional theory (Aziz and Settari 1979, Dullien 1992 as The fractional mobility λ is identical to the harmonic mean mobility λ = λ W λ O /(λ W + λ O ), if the viscous coupling between the percolating phases is negligible (i.e. λ 13 = 0). Inserting equations (28) and (30) into equation (17) gives a set of three coupled nonlinear partial differential equations for the saturations S W , S 2 , S 4 . If the flux Q(t) is imposed externally, it is not necessary to solve the elliptic equation (28) for the pressure, and (31) is a closed system of partial differential equations. The problem is now reduced to finding solutions for equation (31) subject to suitable initial and boundary conditions.

General solution
To avoid overcrowding the water saturation S W with too many indices, the index W is omitted in the subsequent sections. From now on S represents the water saturation instead of S W .

Method of characteristics
Inserting the body forces from equation (10) and the mass exchange terms from equations (5) into (31) one finds where is the gravitational force. The fractional flow function f W and fractional mobility λ are given explicitly by Equations (32b) and (32c) have the solutions (Hilfer 2006c) where are the initial conditions. The limiting saturations S W * , S * 2 , S * 4 are given by equations (6). For imbibition processes (i.e. for ∂ t S > 0) holds for drainage processes (i.e. for ∂ t S < 0). The type of process will at times be indicated by superscripts dr or im. Thus, S im 2 (x, t) is defined by inserting equation (38) into equation (36), and S dr 2 (x, t) is defined by inserting equation (39) into equation (36). Inserting equations (36) with these specifications into (32a), one finds the quasilinear equation where with S 2 = S 2 (S, S W0 , S 20 ) and S 4 = S 4 (S, S W0 , S 40 ) defined by inserting equations (36) into (35). The initial conditions (37) introduce an explicit x-dependence, and the parameter functions Q(t), ϑ(t) an explicit t-dependence. Differentiation of (40) leads to the equation The integral surface S = S(x, t) in (x, t, S)-space solving this quasilinear equation is obtained by solving the autonomous system of ordinary differential equations dt dτ = 1, for the characteristic curves with the curve parameter τ . The functions are considered as functions of S, x and t.

12
The Cauchy problem for (43) consists in finding the solution S(x, t) for given data. The data are written parametrically as where t 0 , x 0 , S 0 are given functions and the map [ξ , , is a parametrization of the Cauchy data. Integrating equations (43) from τ 0 to τ gives with integration constants x 0 (ξ ), S 0 (ξ ) and t 0 (ξ ) given by the problem data.
Of particular interest are data in the (x, S)-plane and their time evolution. They correspond to the special case where t 0 (ξ ) = t 0 for all ξ ∈ [ξ , ξ r ] is a fixed time instant. In this case the choice τ 0 = t 0 leads to τ = t. Eliminating τ the solution at time t > t 0 becomes where ξ ∈ [ξ , ξ r ] is the curve parameter in the (x, S)-plane. This equation suggests that one can interpret g as the velocity of the characteristics and h as the rate of change of the saturation along a characteristic. If the functions x(t, ·) or S(t, ·) can be inverted, i.e. the equations x(t, ξ ) = x or S(t, ξ ) = S can be solved for ξ , at fixed t, then we write by abuse of notation. The saturation profiles S(x, t) defined in equation (48a) may become multivalued. Multivalued saturation profiles are physically not admissible. They indicate the appearance of jump discontinuities in the saturation profile.

Jump discontinuities in the saturations
Location and height of jump discontinuities in the saturation profile S(x, t) are obtained from the inverse function x(S, t) defined in equation (48b). Multivaluedness of S(x, t) corresponds to nonmonotonicity in x(S, t). Nonmonotonicity implies the appearance of local maxima x max and minima x min in x(S, t). The problem of resolving the multivaluedness into jump discontinuities in S(x, t) corresponds to the problem of eliminating the nonmonotonicity in x(S, t) by inserting constant plateaus. The problem is then to determine for every maximum-minimum pair (x (i) such that the conditions x(S (i) are simultaneously fulfilled. For continuous x(S, t) the existence of a solution to this problem is guaranteed by the mean value theorem. Equations (50b)-(50d) are three nonlinear equations for the three unknowns x (i) (t), S (i) In general, the nonlinear equations have to be solved numerically and their solutions are not unique. Once The position where where > 0 are the limiting upstream (+) and downstream (−) flux functions.

Spatial discontinuities in the flux function
A second type of (immobile) discontinuity (suspended shocks) may appear if the initial conditions S W0 (x), S 20 (x), S 40 (x) are discontinuous. Consider a point x = x d in space at which at least one initial saturation S i0 with i ∈ {W, 2, 4} is discontinuous. Now equation (42) does not apply at x = x d and the function h(S, x, t) in equation (44b) is not defined. Mass conservation, expressed in (51), degenerates into flux continuity, because the lefthand side vanishes. Setting the left-hand side in equation (51) to zero requires the condition to hold at all times t > 0. Now the limiting flux functions are defined as with > 0 and the limiting saturations are defined as Equation (53) is an implicit nonlinear condition for the limiting saturations S d − (t) and S d + (t). Depending on the monotonicity of the flux functions it may be fulfilled by more than one pair of limiting saturations. In section 7 spatial discontinuities in monotone flux functions are solved.

Generalized Buckley-Leverett (BL) problem with and without flow reversal
The BL frontal advance theory refers to a hyperbolic approximation, which has become popular, because of its simplicity and versatility in petrophysical applications (Lake 1989). One considers a one-dimensional homogeneous porous medium of infinite extent that is oriented perpendicular to gravity (ϑ = 0 in equation (33)) so that F G = 0 (56) holds.
The initial conditions for the BL problem are specified parametrically as (see equation (47)) where i = W, 2, 4. 3 The boundary conditions at the left and right limits are where i = W, 2, 4. The medium is flooded with the flux protocol where Q 0 > 0 and t r is the instant at which the flow is reversed. Solutions are to be calculated for t ∈ [t 0 , T ] with t 0 < t r < T .

Parameter data
In order to provide a basis for a comparison with experiments, some model parameters have been obtained by fitting to experimental data in Hilfer (2006b). These are listed in table 1. Values for the viscous resistances were not needed in Hilfer (2006b), because only stationary and quasistationary solutions were considered there. Realistic values for viscosities and permeabilities are µ W = µ O = 0.001 kg m −1 s −1 and k = 2.05 × 10 −12 m 2 . Following equations (22) this corresponds to viscous resistances R 11 = R 33 = 5.886 × 10 7 kg m −3 s −1 . The viscous coupling is assumed to occur only between the fluids and the solid, but not between different fluids, so that R 13 = 0, R 14 = 0, R 23 = 0. The method presented here is not restricted to this situation. Including viscous fluid-fluid coupling produces different fronts also ranging from pure rarefactions to pure shocks. Similar to the traditional theory, the fractional flow function f W and the fractional mobility λ in equation (41) are only functions of S. In contrast to the traditional theory however, an infinity of scanning curves is simultaneously available without requiring newly parametrized f W or λ. Figure 1 shows graphs of f W and their derivatives for the four characteristic processes of the hysteresis loop. Primary imbibition is shown as solid curves, secondary imbibition as dashed curves. Primary drainage curves are dash-dotted and secondary drainage curves are dotted. Arrows indicate the direction of the process. The monotonicity of f W is crucial for finding unique solutions in the next section.

Generalized BL frontal advance
In this section the general solution obtained in section 5 with the method of characteristics is applied to solve the initial and boundary value problem from section 6, which is similar to the familiar BL problem (Buckley and Leverett 1942). To do so, we need to discuss spatial disccontinuities in the flux function.

Prerequisite: spatial discontinuities in monotone flux functions
The initial conditions (57) for the BL problem 4 introduce a spatial discontinuity at x = x d = 0 of the type discussed in section 5.3. This requires us to discuss the solution of equation (53) with χ > 0 at both sides of the spatial discontinuity. The limit S 0i + χ is applied for imbibition processes and the limit S 0i − χ for drainage processes. In general the two limits g im i = g dr i differ. Figure 2 illustrates the functions from equations (60). The upper graph shows the constant initial data S , S r for the saturation in the two subdomains separated by a discontinuity. The flux functions G (S), G r (S) are shown in the center of the figure. They are monotone and continuous, but not differentiable. The fluxes at the initial saturations S , S r are marked by circles. The lower graph shows the partial derivatives g (S), g r (S). Note that the two limiting values of g , g r at the initial saturations S and S r are different.
With these preparations the solutions of equation (53) can be obtained. Because gravity is absent, the flux functions are monotonically increasing for t < t r , and hence the slopes g − (S ) > 0, g + (S r ) > 0 are positive. As a result, all characteristics with initial saturation S(x, t = 0) = S arrive at x d from the left (i.e. from the region x < x d ). The characteristics with initial saturation S(x, t = 0) = S r on the other hand also move to the right and cannot reach the spatial discontinuity at x d . It follows that S d − = S for all times. Given S d − = S , the constraint (53) yields an implicit nonlinear equation for S. It has a unique solution, because the flux functions G − (S), G + (S) are monotone in S. Its solution yields the right saturation S d + of the spatial discontinuity as S d + = G −1 + (G − (S )).
4 For simplicity, constant initial saturations S , S 2 , S 4 for x < x d and S r , S 2r , S 4r for x > x d at both sides of the spatial discontinuity are discussed here. Nonconstant initial saturation profiles introduce a time dependence in . This time dependence requires additional considerations but does not introduce conceptual difficulties. Furthermore, it is assumed that the flux functions do not explicitly depend on time. Again, relaxing this constraint requires additional considerations but does not introduce conceptual difficulties.   Open circles show values at S and S r . Open diamonds show the two values at g − (S ) and the plus symbol shows the drainage limit of g dr + (S r ). The values at S d + are illustrated by filled circles. Filled squares show the values at the jump discontinuity. The left inset shows a conceptual picture of the jump discontinuity and the right inset illustrates schematically the drainage and imbibition branch of G − at S . The gray area illustrates constraint (50d) for the jump discontinuity. The dotted horizontal line in the upper figure shows the flux at the spatial discontinuity G d .
The open diamonds on the dashed curve illustrate the values of the two limits g dr − (S ) and g im − (S ). It is not possible to assign a single value to g − (S ), because the saturation at the left side of the spatial discontinuity does not change with time and hence neither imbibition nor drainage takes place there. The filled circle shows (S d + , g + (S d + )). The upstream value g + (S (s) − ) and the downstream value g + (S (s) + ) of the jump discontinuity are depicted as squares and are connected by a dotted line. The equal area construction given in equation (50d) is illustrated by the two gray regions. The region above the dotted line has the same area as the region below it. The open circle on the dash-dotted curve illustrates g im + (S r ), the irrelevant drainage value g dr + (S r ) is shown as a plus symbol. The corresponding saturation profile is illustrated below in figure 4, where also the general solution is derived.

Advancing flow
With these preparations it is now possible to solve the BL initial and boundary value problem formulated in section 6.1. The initial conditions (57) are inserted into equations (41) and (44). In this section the case t < t r is discussed. The case t > t r is addressed in the next section.
For t < t r the flux is constant and positive Q(t) = Q 0 > 0 and G is independent of t where the short forms G (S) = G(S, S , S 2 , S 4 ) and G r (S) = G(S, S r , S 2r , S 4r ) are introduced. Its derivative g with respect to S, defined in (44a), reads g(S, x, t) = g (S), for x < 0, g r (S), for x > 0 (64) in obvious notation. Inserting the piecewise constant initial data (57) into equation (44b) yields At the spatial discontinuity x d = 0 the saturations are related by equation (53). The derivative g(S, x, t) > 0 is positive for all S ∈ [0, 1], for all x ∈ R and for all t ∈ [t 0 , t r ]. Following the previous section 7.1 the saturations at the spatial discontinuity are S d − = S and S d + is determined by equation (62).
For 0 ξ 1 the result is For ξ 1 it reads as If saturation profiles S(x, t) become multivalued, shocks are determined by equation (50). For the problem at hand the number of shocks is zero or one. A shock is denoted (i) = (s) when occurring. The shock occurs in the domain x > 0, because in the region x < 0 the saturation is constant for all time. The downstream saturation value of the shock S (s) + = S r and the upstream saturation S (s) − is independent of time, because of the initial conditions (57). Inserting this into equation (67a) the position of the shock becomes x (s) (t) = tv (s) with v (s) = g r (S (s) − ).

Inserting this and equations (66)-(68) into (50) one obtains an implicit equation for the upstream saturation S (s)
− at the jump discontinuity Dividing equations (66a), (67a) and (68a) by t it becomes apparent that the saturation S(x, t) is self-similar in x/t and reads (70) Figure 4 shows a conceptual picture of the solution at an instant t. The saturation profile is depicted as a thick solid line. The spatial discontinuity in the initial conditions is indicated by a vertical line. Thin dotted curves conceptually show g −1 and g −1 r . Characteristic values are illustrated by the same symbols that were used in figure 3. The figure illustrates how the constant saturation at the left is connected via a jump at the spatial discontinuity and a rarefaction shock with the constant saturation at the right.
Several displacement processes are discussed below. Table 2 lists the initial conditions used for the displacement processes which are illustrated in figures 5 and 6. Characteristic values of the solutions to these problems are summarized in table 3. The parameters are listed in table 1 above. Figure 5(a) shows a primary imbibition process. A primary imbibition is a displacement of the nonwetting fluid by the wetting fluid in a completely dry porous medium. The initial saturations for the primary process hence are S = 1 and S r = 0. The initial saturations of the nonpercolating phases are S 2 = 0, S 2r = 0, S 4 = 0, and S 4r = 0. The figure shows the initial saturation jump at x = 0 as a dashed line. After t = 15d a water shock front has invaded the porous medium and has reached x (s) (15d) = 22.23 m. It is followed by a rarefaction fan spanning the interval between x = 0 and the shock front. During the displacement the oil phase Table 2. Initial values for displacement processes illustrated in figures 5 and 6.     Figure (a) shows profiles for a primary imbibition process and figure (b) shows profiles for a secondary imbibition process. Profiles are shown at t = 15d as solid curves. The initial conditions are listed in table 2 and depicted as dashed lines. Note, that some curves coincide with the x-axis.
may break up and nonpercolating oil remains immobile in the medium. The water phase remains continuous, because it is invading. Figure 5(b) shows a secondary imbibition process. A secondary imbibition is a displacement of the nonwetting fluid by the wetting fluid in a medium containing only irreducible water at the beginning. Hence the initial conditions are S = 1 − S O im , S r = S W dr , S 2 = 0, S 2r = S W dr , S 4 = S O im , and S 4r = 0. The figure shows the initial conditions as dashed lines. After t = 15d a water shock front has invaded the porous medium and has reached x (s) (15d) = 30.07 m. It is followed by a rarefaction fan. During the displacement the oil phase eventually breaks up and nonpercolating oil remains trapped inside the medium. The nonpercolating water phase on the other hand coalesces with the invading percolating water phase.
A comparison of both results shows that the shock of the secondary imbibition process is significantly faster than the one of the primary process whereas the shock height is a little bit smaller. See table 3 for numerical values. The higher velocity is due to the coalescence of initially present nonpercolating water with the imbibing percolating water. The coalescence increases the mobility λ W and the flux Q W of the water. Figure 6 illustrates the influence of the mass exchange efficiency parameters η 2 , η 4 on the saturation profiles. A comparison of profiles for different η 2 , η 4 in a secondary imbibition process is shown. Solid curves represent saturation profiles for η 2 = 4, η 4 = 3, dashed curves for η 2 = 1, η 4 = 1 and dotted curves for η 2 = 10, η 4 = 10. The other parameters are chosen identical to the preceding examples. The parameters η 2 , η 4 significantly change shock velocities v (s) and shock heights S (s) − . For larger η 2 , η 4 the shock velocities are larger. The shock heights are correspondingly smaller, because the total fluid volume is conserved. These solutions have also been checked numerically using an adaptive moving grid algorithm (Doster 2011). Figures 7(a) and (b) show two different displacement processes in a homogeneous porous medium with constant initial water saturations. The subdomains x < 0 and x > 0 differ only in the ratio between percolating and nonpercolating phases, but their total water saturation S W is identical. Initially S W = 0.6 in the whole domain. In figure 7(a) the region x < 0 has been prepared through a primary drainage, while x > 0 has been prepared through a primary (b) (a) Figure 7. Saturation profiles for a BL-type problem with constant initial water saturation S = 0.6 in the whole domain. Both subdomains x < 0 and x > 0 differ only in the ratio between percolating and nonpercolating phases. Figure (a) shows profiles for initial conditions where the left domain has been prepared through a primary drainage and the right domain through a primary imbibition. Figure (b) shows profiles for mirror reflected initial conditions, i.e. the left domain is prepared through a primary imbibition and the right domain through a primary drainage. Numerical values are listed in table 2. Initial conditions are depicted as dashed lines and profiles at t = 15d as solid lines.
imbibition. In figure 7(b) the initial conditions are mirror reflected. The initial conditions are depicted as dashed lines and solutions at t = 15d as solid lines.
When water is injected from the left one finds S d + < S r in figure 7(a). Thus, a drainage shock is moving from the left to the right, reaching x (s) (15d) = 60.02 m after t = 15d. Drainage occurs, because on the left side (x < 0) the amount of mobile, i.e. percolating, water is less than on the right side. The water flux Q W on the left is less than the flux at the right side, and the water flux which is required for a constant water saturation in the right side is not provided by the left side. The opposite holds true for figure 7(b). Here the higher water flux from the left side leads to an increase of the water saturation at the right S d + > S r and an imbibition wave is propagating into the right domain. It is a pure rarefaction. Both examples illustrated in figure 7 differ from and go beyond the traditional theory, where the saturation would remain unchanged at its initial value S W = 0.6 throughout the whole domain.

Reversed flow
We emphasize that our extended theory predicts hysteresis. This will be illustrated here by reversing the flow, i.e.
at a time instant t r . Hence, the front will move back and the character of the process changes locally. An imbibition now becomes a drainage and a drainage becomes an imbibition. Contrary to the traditional theory, the full hysteresis loop and all scanning curves are automatically contained in our extended theory. An imbibition followed by a drainage will be discussed below without loss of generality. For drainage followed by imbibition only the wording has to be changed accordingly.
The solution profiles S(x, t r ), S 2 (x, t r ), S 4 (x, t r ) given in equation (70) are the initial conditions for the reversed flow problem. For t t r the domain is generally divided into four domains. The left domain x < 0 is denoted by an index . The second domain 0 < x < t r g r (S d + ) is denoted by an index c, because the saturation at t = t r is constant there. The third domain t r g r (S d + ) < x < x (s) (t r ) corresponds to the rarefaction of solution (70) at t = t r . It will be denoted by an index m. The right domain x > x (s) (t r ) will be denoted by an index r . Assuming t r = 15d and the initial conditions 6 from table 2 the solid curves in figure 6 show the initial conditions for a reversed flow. Then the domains are : x ∈ [−∞, 0]m, c : x ∈ [0, 11.68]m, m : x ∈ [11.68, 34.08]m and r : x ∈ [34.08, ∞]m.
The saturations in , c and r are constant at t = t r . Hence, the flux functions G , G c , G r do not explicitly depend on x and sources in the quasilinear partial differential equation (42) are absent so that holds. The current flows from right to left and thus the saturation is constant for all t > t r and x > x (s) (t r ). At x = x (s) (t r ) there is a discontinuity in the initial conditions because of the shock. This leads to a second discontinuity in the flux function which is denoted by a superscript ds. The constraint (53) holds at the spatial discontinuity and g(S, x) < 0 for all S ∈ [0, 1], for all x ∈ R and for all t ∈ [t r , T ]. Following section 7.1 the saturations at the spatial discontinuity are S ds + = S r and S − is the solution of equation (62) with S d + = S r . The interface at x = t r g r (S d + ) does not require special treatment and the spatial discontinuity at x = 0 in the initial conditions is due to the spatial discontinuity in the flux functions during the frontal advancement for t < t r and is treated analogously. In the region called m a source term is appearing in the quasilinear partial differential equation (42) h (S, x, t) for t r g r (S d + ) < x < x (s) (t r ) and t > t r . The solution of the advancing flow (70) is now parametrized with the parameter ξ ∈ R as and to specify the Cauchy data at time t = t r . Inserting this into equation (47) again yields solutions for the different parameter domains. In contrast to the previously discussed Riemann problem, the initial conditions for the drainage are not homogeneous and an analytical solution is not obvious. However, the set of ordinary differential equations (43) can be solved numerically. We chose an adaptive Runge-Kutta method of fifth respectively fourth order provided within the software package Matlab as ode45 (Dormand and Prince 1980).
The solution in (x, t, S)-space may yield a multivalued function S(x, t). The multivaluedness is resolved by equations (50). Here equations (50) are solved numerically by a bisection method.
As an example, a secondary imbibition followed by a drainage is illustrated in figure 8. Note, that in this case S d + = S and g r (S d + ) = 0. The domain c shrinks to a point x = x (s) (t r ). Figure 8(a) shows profiles at equidistant time steps for t t r . The dashed profiles are the initial conditions. Figure 8 (b) shows profiles at equidistant time steps for t t r . The dashed profiles correspond to t = t r . The time step between adjacent profiles in both subfigures is t = 3d = t r /5. For 0 t t r the current flows from left to right. For x < 0 m the saturation is constant and for x > 0 m a secondary imbibition process takes place. For t t r the current is reversed and flows from right to left. A rarefaction shock then propagates from right to left. Contrary to before the saturation S (s) + (t) ahead of the shock is not spatially constant. The upstream saturation S (s) − (t) behind the shock and the velocity v (s) (t) are also changing with time. When the shock returns to x = 0.0m it turns into a standard rarefaction shock, whose shock velocity, up-and downstream saturations are the constant values of a secondary drainage. If one wanted to treat this problem with the traditional theory, this would require an infinite number of different scanning curves (in the domain m). This is not the case for our theory.

Conclusion
The present paper contributes to the long-standing open problem of formulating a macroscopic hydrodynamic theory of capillary hysteresis and fluid entrapment in porous media. Wider implications of our results beyond physical fluid dynamics and wetting in porous media concern nonequilibrium physics, where hysteresis and memory effects are ubiquitous, but poorly understood phenomena. Mathematically, our work provides analytical and numerical solutions for systems of nonlinear partial differential equations. Moreover, our results are important for many applications, such as enhanced oil recovery, CO 2 -sequestration, fuel cells, soil physics, or hydrology.
In summary, this paper has advanced the traditional, 70-year-old, widely used BL theory for two-phase flow in porous media by taking into account the differences of percolating and nonpercolating fluid parts. The mathematical model is based on extending the state space by two additional saturations and velocities. A fractional flow formulation of the extended theory was given. A hyperbolic limit of the extended theory can be formulated by assuming immobile nonpercolating phases and equal pressures for the percolating fluid phases. The method of characteristics was applied in this limit to solve the BL problem. Because of hysteresis, the initial and boundary value problem gives rise to spatially discontinuous flux functions, even for homogeneous media.
The solutions of the generalized BL problem show significant differences of shock heights and speeds between primary and secondary processes. It was found that mass exchange parameters influence the velocities and heights of the shocks. It was further shown that even though the initial water saturation is constant, drainage and imbibition waves may form due to differences in the hysteresis loop. The ability to switch dynamically between drainage and imbibition processes within our extended theory permits solutions also for reversed flow. The extended theory does not require different parameter sets to treat all these cases. Drainage and imbibition curves do not have to be selected in advance. This is essential for many experimental and technological displacement processes, where it is impossible to know the processes from the outset.
For the benefit and convenience of the general reader we briefly recall traditional BL theory. The material in this appendix is not new (see e.g. Buckley and Leverett (1942), Aziz and Settari (1979), Chen et al (2006), Hilfer (2006bHilfer ( , 2006c). We emphasize the close analogy between traditional BL theory and our generalization.

A.1. Macroscopic phase structure
The pore space is called P, the matrix is called M, and the sample is S = P ∪ M. The traditional macroscopic theory distinguishes two phases: a wetting phase, called water and denoted W, and a nonwetting phase, called oil and denoted O. The microscopic fluid configuration (on the pore scale) is not resolved, and the phases W and O are viewed as being simultaneously present at each macroscopic position x.

A.2. Balance laws for mass, momentum and volume
The traditional theory starts from the fundamental balance laws of continuum mechanics for water W and oil O. The law of mass balance in differential form reads (cf. equation (2)) , v i (x, t) denote mass density, volume fraction and velocity of phase i = W, O as functions of position x ∈ S ⊂ R 3 and time t ∈ R + . The mass transfer rates M i give the amount of mass by which phase i changes per unit time and volume. They can be used to model mass sources or chemical reactions. Momentum balance for the two fluids requires in addition (analogously to equation (7)) where i is the stress tensor in the ith phase. F i is the body force per unit volume acting on the ith phase, m i is the momentum transfer into phase i from all the other phases, and D i /Dt = ∂/∂t + v i · ∇ denotes the material derivative for phase i = W, O.
Defining the saturations S i (x, t) as the volume fraction of pore space P filled with phase i one has the relation φ i = φ S i where φ is the porosity of the sample. Expressing volume conservation φ W + φ O = φ in terms of saturations yields (analogously to equation (1)) In order to get the traditional theory these balance laws for mass, momentum and volume have to be combined with specific constitutive assumptions for M i , m i , F i and i .

A.3. General constitutive assumptions
As a first approximation it is usually assumed that the porous medium has a solid, ideally rigid matrix that is also macroscopically homogeneous (analogous to equation (3)) φ(x, t) = φ = const (A.4) although this assumption can be relaxed, and rarely holds in practice (Hilfer and Helmig 2004). Let us further assume that the fluids are incompressible so that (analogous to equation (4)) where the constants W , O are independent of x and t. One assumes next that the stress tensor of the fluids is diagonal (see equations (9a) and (9c)) where P W , P O are the fluid pressures. Realistic subsurface flows have low Reynolds numbers so that (analogously to equation (8)) the inertial term can be neglected in the momentum balance equation (7). It is further assumed that the body forces are given by gravity (as in equation (10) with ϑ = π/2). As long as there are no sources or sinks and no chemical reactions between the fluids, the mass transfer rates vanish, so that holds.

A.4. Viscous drag
Momentum transfer between the fluids and the rigid walls of the porous medium is assumed to be governed by viscous drag in the form (analogously to equation (11)) where µ W , µ O are the constant fluid viscosities, k is the absolute permeability, and k r W (S W ), k r O (S W ) are the so-called relative permeabilities functions of water and oil. Note that these functions are in general nonlinear, and partly responsible for the nonlinearities of the theory.

A.5. Capillarity
Inserting the constitutive assumptions into the general balance laws (A.1)-(A.3) yields nine equations for ten unknowns S W , S O , P W P O , v W , v O . An additional equation is needed. Based on observations of capillary rise in regular packings (Smith et al 1931) it was argued in Leverett (1941) that the pressure difference between oil and water should, in general, depend only upon saturation where σ WO is the oil-water interfacial tension and κ(S W ) is the mean curvature of the oil-water interface. This assumption has remained the cornerstone of the theory of macroscopic capillarity for 60 years, and it is being challenged here. The nonlinear constitutive parameter function P c (S W ) is called the capillary pressure-saturation relation and it is supposed to describe the macroscopic effect of capillarity in hydrostatic equilibrium (without flow).

A.6. Ad hoc assumptions
The functions P c (S W ) and k r W (S W ), k r O (S W ) are treated as constitutive functions that are supposed to represent 'rock properties' (see e.g. chapter 3 in Chen et al (2006)). It is well known, however, that these functions are not intrinsic constitutive parameters, but depend on time, position, fluid configurations, and on the type of flow process. Technically, they are fit functions that are chosen freely to approximate direct or inverse experiments. Fit functions can often be subsumed under the general form and P † c , P ‡ c , k r W † , k r O † , a W , a O , S Wi , S Or and ν 1 , . . . , ν 11 are 19 parameters. Different forms using logarithms or exponential functions have also been proposed. The parameters depend on the type of displacement process. Frequently the number of parameters is reduced on the basis of observations as well as calculations for capillary tube models. A popular choice for primary drainage assumes S Or = 0, P ‡ c = 0, a W = 0, a O = 0, ν 5 = 0, ν 9 = 0 leaving seven parameters free.

A.7. BL approximation
The BL approximation firstly neglects gravity Secondly, it ignores capillarity effects by assuming that (analogous to equation (25)) .15) and hence P c (S W ) = 0 holds 5 . This assumption is justified by practitioners in petroleum engineering by the argument that capillary phenomena are negligible on large scales of hectometers or kilometers. Combining all these assumptions allows one to formulate the theory as an elliptic equation for the pressure plus the so-called BL equation (Buckley and Leverett 1942) (in one dimension) for the water saturation, where is the BL fractional flow function. The total volume flux (analogous to equation (27)) is independent of position and assumed to be given as a function of time (usually constant). The BL equation is a nonlinear partial differential equations. Smooth initial data may evolve into a shock front.