A vanishing dynamic capillarity limit equation with discontinuous flux

We prove existence and uniqueness of a solution to the Cauchy problem corresponding to the dynamics capillarity equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\left\{ \begin{array}{ll} \partial _t u_{\varepsilon ,\delta } +\mathrm {div} {\mathfrak f}_{\varepsilon ,\delta }(\mathbf{x}, u_{\varepsilon ,\delta })=\varepsilon \Delta u_{\varepsilon ,\delta }+\delta (\varepsilon ) \partial _t \Delta u_{\varepsilon ,\delta }, \ \ \mathbf{x} \in M, \ \ t\ge 0\\ u|_{t=0}=u_0(\mathbf{x}). \end{array}\right. } \end{aligned}$$\end{document}∂tuε,δ+divfε,δ(x,uε,δ)=εΔuε,δ+δ(ε)∂tΔuε,δ,x∈M,t≥0u|t=0=u0(x).Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathfrak {f}}}_{\varepsilon ,\delta }$$\end{document}fε,δ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_0$$\end{document}u0 are smooth functions while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}ε and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta =\delta (\varepsilon )$$\end{document}δ=δ(ε) are fixed constants. Assuming \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathfrak {f}}}_{\varepsilon ,\delta } \rightarrow {{\mathfrak {f}}}\in L^p( {\mathbb {R}}^d\times {\mathbb {R}};{\mathbb {R}}^d)$$\end{document}fε,δ→f∈Lp(Rd×R;Rd) for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1<p<\infty $$\end{document}1<p<∞, strongly as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \rightarrow 0$$\end{document}ε→0, we prove that, under an appropriate relationship between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}ε and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta (\varepsilon )$$\end{document}δ(ε) depending on the regularity of the flux \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathfrak {f}}}$$\end{document}f, the sequence of solutions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(u_{\varepsilon ,\delta })$$\end{document}(uε,δ) strongly converges in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^1_{loc}({\mathbb {R}}^+\times {\mathbb {R}}^d)$$\end{document}Lloc1(R+×Rd) toward a solution to the conservation law \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \partial _t u +\mathrm {div} {{\mathfrak {f}}}(\mathbf{x}, u)=0. \end{aligned}$$\end{document}∂tu+divf(x,u)=0.The main tools employed in the proof are the Leray–Schauder fixed point theorem for the first part and reduction to the kinetic formulation combined with recent results in the velocity averaging theory for the second. These results have the potential to generate a stable semigroup of solutions to the underlying scalar conservation laws different from the Kruzhkov entropy solutions concept.


Introduction and notation
Flow in a two-phase porous medium is governed by the Darcy law [7] q = −K(S) (∇p + ρge d ) , (1.1) where e d = (0, . . . , 0, 1), is the direction of gravity. The quantity S is the saturation, p is the pressure, and (the vector) q is the flow velocity of the wetting phase (usually water, while the non-wetting one is oil or a gas). The Darcy law represents conservation of momentum and, in order to close the system, we also need the conservation of mass ∂ t S + divq = 0. (1.2) is monotonic (see, e.g., the introduction of [53]). Therefore, many attempts to explain the gap between the standard theory (provided by the Buckley-Leverett or Richards [45] equations) have been put forth recently. A purely mathematical approach can be found in [27], where δ-type solutions to (1.3) are constructed and δ-distributions appearing as a part of the solution are explained as an inadequacy of the model. In [13], a model is suggested that involves the notion of non-local energy, which eventually leads to a fourth-order equation whose solution (with the Riemann initial data given above) has a shape corresponding to the experimental results from [14].
In [53] one can find an interesting approach explaining the discrepancy between theory and experiment described above, based on the dynamics capillarity concept introduced in [22,23], which has drawn a lot of interest (especially after the publication of [53]; see also [1]). Namely, in [22,23] it was supposed that the capillary pressure depends not only on the saturation S, but also on the time derivative of the saturation (of the wetting phase): Taking this into account and proceeding along the lines of deriving the Buckley-Leverett equation, one reaches a nonlinear pseudo-parabolic equation ( [53, (1.16)]) which, after linearization of higher-order terms reduces to The x-dependence of the flux means that we assume that the medium in which we consider the phenomenon is heterogeneous, i.e., that it has different properties at different points (e.g., since in some parts of the medium we have sand while in some others clay). Moreover, we shall assume that the flux is discontinuous with respect to the space variable x, which means that the medium experiences abrupt (discontinuous) changes in its properties (one can imagine that we have sand which is highly permeable adjacent to the clay layer, which is weakly permeable-permeability is discontinuous in such a medium). Let us remark here that, stipulated by different applications, evolutionary equations with discontinuous flux have attracted considerable attention recently. For a (non-exhaustive) selection of recent results, cf. [2,3,12,29,43] and references therein.
Explicitly, we consider the conservation law ∂ t u + divf(x, u) = 0, (1.6) where and ∂ λ f i ∈ L p (R d × R) for some 1 < p < ∞ and all i = 1, . . . , d; is the space of Radon measures, and there exists a constant β > 0 and a finite Radon measure μ ∈ M(R d ) such that (in the sense of measures) We then perturb (1.6) by (a linearized) vanishing dynamic capillarity limit: is a bounded family of compactly supported functions equal to one on the ball B(0, 1/ε) ⊂ R d+1 and such that ∇K ε L ∞ (R d+1 ) ≤ 1. We remark here that (1.7) reminds on the diffusion-dispersion regularization (see, e.g., [4,31]), but it contains t-derivative in the third-order term which makes significant difference between the two situations. We will comment on this in more details later. More precisely, we suppose that (1) , ω (2) test functions with unit integral on R d and R, respectively, and with an appropriate n(ε) such that n(ε) → 0 as ε → 0 (the form of n(ε) will be precisely determined later). Note that We supplement (1.6) with the initial data while we take where, for every i = 1, . . . , d, n(ε) is as above, satisfies (1.11). We shall, however, use (1.11) as an assumption.
With regard to pseudo-parabolic equations, we have already explained their importance in the porous media theory. Besides a possible justification of the experimental results from [14], there is a confirmed application to the seepage of homogeneous fluids through a fissured rock [5]. Also, they describe the unidirectional propagation of nonlinear, dispersive long waves [6,51] (where u is typically the amplitude or velocity), or population dynamics [41] (where u represents the population density). In [48,52], the authors investigated the initial-boundary value problem and the Cauchy problem for linear pseudo-parabolic equations and established the existence and uniqueness of solutions. As for the nonlinear variants, one can find numerous results even for singular pseudo-parabolic equations and degenerate pseudo-parabolic equations (see, e.g., [9,28,30,34,39,44] and the references therein). Together with the existence and uniqueness results, among the given references, one can also find properties of solutions, such as asymptotic behavior and regularity. However, we are not aware of any corresponding results for the Cauchy problem for an equation of type (1.5).
This problem will be considered in Sect. 2. We shall use an approach that is characteristic for the theory of wave equations [47]: we shall apply the Fourier transform with respect to x, solve the ordinary differential equation so obtained with respect to t and prove the following theorem (the constants ε and δ are omitted for simplicity): There exists a unique solution to In the next step (Sect. 3), we are going to let the perturbation parameter ε → 0 in (1.7). This type of problem for conservation laws (when there is not only vanishing viscosity, but also third-or higher-order perturbation) was first addressed in [46] and received considerable attention after that. A thorough analysis of this kind of limit in the one-dimensional situation with a regular flux can be found in [38], where the theory of non-classical shocks for conservation laws was essentially initiated. The standard approach here is to rewrite the equation under consideration in the kinetic formulation or using Young measures (which is essentially equivalent), and then applying velocity averaging results [4,26], compensated compactness (in one-dimensional situations) [10,11,46], or Di Perna techniques involving Young measures [15], as done in [31] and many others. Again we note that this list of citations is far from complete. The problem in our case is the low regularity of the flux function f (it can be discontinuous with respect to x ∈ R d ) as well as the multidimensional character of our problem, which prevents us from using the Young measures and compensated compactness approach here (see item (iii) below and Remark 3.9 for the case of a regular flux). We also remark that, to the best of our knowledge, the only diffusion-dispersion type result for equations with discontinuous coefficients is [24].
Moreover, unlike the situation that is typical in the case of the diffusion-dispersion limit ( [4,31,46] etc.), where the existence of a solution to the perturbed equation is assumed together with all necessary properties of the solution, we have proved in the previous section the existence of the functional sequence whose convergence we analyze.
As for the velocity averaging theory, most of the results are given in the case of a homogeneous flux [16,21,42,50] or a flux for which p ≥ 2 [19,35]. We have recently proved [36] the velocity averaging lemma in the case p > 1 and this will enable us to prove the strong convergence of the sequence (u ε ) of solutions to (1.7), (1.10) along a subsequence toward a weak solution to (1.6). The following statements are the main results of the paper: (1.14) as ε → 0, then the family (u ε ) contains a strongly converging subsequence in L 1 loc (R + × R d ) under a suitable non-degeneracy condition (see Definition 3.2); (ii) If δ = O(ε 2 ) and √ δ √ εn(ε) d/2+2 → 0 as ε → 0, and the λ-derivative of the flux ∂ λ f is bounded in addition to condition (C3), then the family (u ε ) is strongly precompact in L 1 loc (R + × R d ) under the non-degeneracy condition; , then the family (u ε ) converges strongly in L 1 loc (R + × R d ) toward the entropy solution to (1.6), (1.9) (no non-degeneracy condition is needed).
We conclude the introduction by noting that the subsequence of (u ε ) given in (i) does not converge toward an entropy admissible solution to the underlying conservation law since its flux is not smooth and for such fluxes the well-posedness theory is not developed yet in full generality (compare with item (iii) above).
Also, in the case (ii), even when the flux is smooth, we do not necessarily have convergence toward the Kruzhkov entropy solution (see [53] for stepwise initial data). On the other hand, the approximating procedure actually models a physical situation in which diffusion and capillary limits have equal effect on the process and, as any physical phenomenon in the macro-world is well-posed, it is therefore expected to generate a stable semigroup of solutions to the underlying conservation law. We will deal with this question in a future research.
The paper is organized as follows: In Sect. 2, we prove existence and uniqueness for the pseudoparabolic problem (1.12), (1.13). In Sect. 3, we derive necessary estimates from (1.7), (1.10) and use it to prove the strong convergence toward a weak solution to the underlying conservation law.

Existence and uniqueness of the solution to the pseudo-parabolic equation (1.12) with (1.13)
Throughout this section, we suppose that f is smooth and compactly supported. Under this assumption, we want to show that (1.12) with the initial condition (1.13) has a unique solution. The strategy for solving this problem is to define a mapping T : represents a solution to with the initial conditions (1.13). Generally, for u 0 ∈ L 2 (R d ), we are seeking weak solutions u of this initial value problem in the sense that, for any smooth test function ϕ with compact support in [0, T ), we require Then, we shall prove that the mapping T possesses a fixed point, which will turn out to be the solution to (1.12), (1.13).
To this end, we need the following consequence of the Leray-Schauder fixed point theorem (cf. [20, Th. 11.3]): Theorem 2.1. Let T be a compact mapping of a Banach space B into itself and suppose that there exists a constant C such that for all u ∈ B and σ ∈ [0, 1] satisfying u = σT u. Then, T has a fixed point, that is, T u = u for some u ∈ B.
Finally, let us fix our conventions for Fourier transform and inverse Fourier transform.
and the inverse Fourier transform is given by where · denotes the scalar product and i is the imaginary unit.
We begin by proving existence and uniqueness of solutions to (2.2), (1.13).
Proof. After applying the Fourier transform with respect to x ∈ R d to (2.2), we obtain We rewrite the last equation in the form It follows from this that any weak solution to (2.2) has t-derivative in L 2 as well, hence in particular is continuous with respect to t and therefore possesses a classical trace on t = 0. Solving the ODE with the initial dataû(0, ξ) =û 0 (ξ), we arrive at . By finding the inverse Fourier transform here with respect to ξ ∈ R d , we indeed obtain a weak solution, so combined with the above considerations both existence and uniqueness follow.
Proof. Let u be the solution to ( where in the last step we used the Plancherel theorem again. From here, since f is compactly supported with respect to x ∈ R d and bounded, we conclude that As for the t-derivative, we have and therefore, repeating the procedure giving us the estimates on ∂ xj u, we get Based on this, we can proceed to applying Theorem 2.1 to prove existence of a solution to (1.12), (1.13).

Theorem 2.4.
There exists a unique solution to (1.12), (1.13). This solution belongs to Proof. Let us first take a sequence of balls B n = B(0, n) centered at 0 and of radius n. Set u n where T is the operator from (2.1) with the initial data u 0 χ Bn and we extend vχ Bn by zero so it is defined on [0, T ] × R d . This map is continuous since T is continuous. In fact, by the proof of Theorem 2.3, T is in fact continuous as a map from . So, we can use Theorem 2.3 and the Rellich theorem (keeping in mind that B n is a bounded set) to conclude that T n is a compact mapping. If we apply Fourier transform with respect to x ∈ R d here and solve the ODE so obtained, we arrive at an expression analogous to (2.4): and repeating the procedure from (2.5) and (2.6), we conclude that (2.3) holds for the mapping T n (keeping in mind that n is fixed).
Thus, for every n ∈ N there exists a function u n ∈ L 2 ([0, T ] × B n ) solving Extending u n by zero outside of B n , we thereby obtain a sequence (u n ) that is bounded in L 2 ([0, T ]×R d ), and locally bounded in H 1 ((0, T ) × R d ). Thus, (u n ) admits a subsequence that locally converges toward It is clear that the function u is the weak solution to (1.12), (1.13) and thus it must belong at least to Moreover, the proof of Theorem 2.3 with v = u is a classical example of a bootstrapping procedure, which enables us to conclude that u actually belongs to each H k ((0, T ) × R d ) (k ∈ N), hence is smooth (due to smoothness of u 0 and f).
Now we turn to the proof of uniqueness. Assume that u 1 and u 2 are solutions to (1.12), supplemented with the initial conditions u 1 t=0 = u 10 and u 2 t=0 = u 20 . After applying the procedure from Theorem 2.3 we conclude that wheref Taking the L 2 (R d )-norm of (2.7), we get, due to (2.8): Applying Gronwall's inequality, we arrive at From here, we see that u 10 = u 20 entails u 1 = u 2 . This concludes the proof.

Vanishing capillarity limit
In this section, we inspect the vanishing capillarity limit of (1.7). By Theorem 1.1, under the assumptions (C1)-(C3) from Sect. 1, the Eq. (1.7) with initial data (1.10) (satisfying (1.11)) possesses a unique solution As announced in Sect. 1, our main result then is as follows:  Below we give a complete proof of (i). Although non-trivial, the other two claims are more standard and we shall only comment on the necessary modifications of the argument for establishing (i).

Definition 3.2. The flux appearing in (1.6) is called non-degenerate if for the function
is not zero on any set of positive measure.
Then, for any ρ ∈ L 2 c (R), there exists a subsequence (h n k ) of (h n ) such that as k → ∞. Now, we proceed to the estimates sufficient to apply the mentioned functional analytic tools. The solution u ε of (1.7) satisfies the following a priori estimates:

Proof. Let us consider a family (u ε ) of smooth solutions to (1.7). Multiplying the equation by
where q ε is defined by (recall that f ε = (f ε 1 , . . . , f ε d )) and normalized by the condition q ε j (x, 0) = 0, j = 1, . . . , d. Integrating over x ∈ R d and 0 ≤ t ≤ T and using integration by parts, we get The relations (3.4), (3.5) and (3.6) are easy consequences of (3.10). To show (3.7), let us now differentiate equation (1.7) with respect to x i and then multiply it by 2∂ xi u ε . Since the right-hand side of (1.7) is linear in u ε , the calculation for that side will remain the same and after integrating the equation over [0, t] × R d and applying integration by parts, we get . Multiplying the above inequality by δ, we have , where in the last inequality we used Eq. (3.10). From here, (3.7) immediately follows.
We next aim to prove inequality (3.8). To this end, we multiply equation (1.7) by ∂ t u ε and then proceed similarly to the above. Namely, and we integrate equation (3.12) with respect to time and space. Integration by part gives i.e., (3.13) By multiplying (3.13) by δ and taking into account (3.5), we conclude which gives the inequality (3.8).
We shall now derive the kinetic formulation for (1.7) and then show that we can apply the velocity averaging result on this form of the equation.  (Kinetic formulation of (1.7)) Assume that the function u ε satisfies (1.7). Then, the function h ε (t, x, λ) = sgn(u ε (t, x) − λ) solves a linear PDE of the form (ii) (G ε 2 ) is bounded in the space of measures M loc (R + × R d × R) and thus strongly precompact in Proof. Before we start, notice that we require merely local convergence (see conditions (i) and (ii) of the lemma), which means that we can always choose ε small enough so that K ε ≡ 1 on the fixed relatively compact set on which we derive the estimates. Therefore, it is enough to prove the lemma under the simplifying assumption Also, in order to avoid proliferation of symbols, we will notationally adhere to integration over the entire space R d × R, with the implicit understanding we are in fact integrating over the support of the corresponding test function (see (3.19)). Then, note that for any To proceed, we multiply equation (1.7) by η (u ε ). We get Now, we apply (3.16) and (3.17) and rewrite the latter equation in the form If we rewrite the latter in the variational formulation, we get for every test function φ ∈ C 2 we arrive at (3.15). Next, we prove (i) and (ii). To begin with, we show that divΓ ε We have for any φ ∈ C 1 c (R + × R d × R) supported in the hypercube with side length M centered at zero: (3.21) To proceed further, we use (3.5). For estimating the term div(f ε ) L 1 (R d+1 ) appearing in (3.5), we use approximations of the form f ε := f ω n(ε) as in (1.8) with n(ε) satisfying (1.14). Since δ = o(ε 2 ), the f ε will converge to f as ε → 0 in L r loc (R d × R) for any 1 r p. Using (C2) and (C3), we get (3.22) where here and below, C denotes a generic constant that may alter from line to line. Thus (keeping in mind (3.5) and (1.11)), from (3.22) we conclude by (1.14).
According to the previous theorem, we see that we can apply the velocity averaging lemma Theorem 3.3. Indeed, we have: Lemma 3.6. Any solution h ε of (3.15) satisfies )). On the other hand, according to condition (C2), we see that (h ε divf ε (x, λ))) is bounded in the space of measures and thus strongly precompact in W −1,r loc (R + × R d × R). This concludes the proof. Now, we are ready to use Theorem 3.3. In fact, we may choose r ∈ 1, d+2 d−1 and s 2 in such a way that and then setp := (1 + 1 s − 1 r ) −1 . Then, fixingp such that 1 p + 1 p = 1, the assumptions of Theorem 3.1 precisely allow to apply Theorem 3.3 (noting that (h ε ) is bounded in any L s loc , hence contains a weakly convergent subsequence). More precisely, one easily checks that anyp ∈ 2d+4 d+4 , ∞ , and thereby anȳ p ∈ 1, 2 + 4 d can be obtained in this way. This is the reason for the specific assumption in Theorem 3.1. Consequently, for any sequence ε n 0 and any ρ ∈ L 2 (R), setting h n := h εn and u n := u εn , the sequence R ρ(λ)h n (t, x, λ)dλ is strongly precompact in L 1 loc (R + × R d ). As we shall see in the next theorem, this implies strong convergence of the sequence (u n ). We first need the following auxiliary result. With the above notations, we finally arrive at: Theorem 3.8. Under the non-degeneracy conditions (3.1), the sequence (u n = u εn ) of solutions to (1.7), (1.10) strongly converges along a subsequence in L 1 loc (R + × R d ).
We actually proved only item (i) of Theorem 3.1. The other two items (item (ii) and item (iii)) can be proven by an adaptation of the proof for item (i). We provide a more precise explanation in the following remark. Remark 3.9. Derivation of (ii) and (iii) from Theorem 3.1.
If we additionally assume that ∂ λ f ≤ C < ∞, then we can use (3.8) in the case δ = O(ε 2 ) and √ δ √ εn(ε) d/2+2 → 0 and the considerations thereafter will remain the same. This gives (ii). Indeed, we have to check whether the sequence of equations (3.30) satisfies the conditions of Theorem 3.3. To this end,