Piecewise Constant Subsolutions for the Muskat Problem

We show the existence of infinitely many admissible weak solutions for the incompressible porous media equations for all Muskat-type initial data with C3,α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C^{3,\alpha}}$$\end{document}-regularity of the interface in the unstable regime and for all non-horizontal data with C3,α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C^{3,\alpha}}$$\end{document}-regularity in the stable regime. Our approach involves constructing admissible subsolutions with piecewise constant densities. This allows us to give a rather short proof where it suffices to calculate the velocity and acceleration at time zero - thus emphasizing the instantaneous nature of non-uniqueness due to discontinuities in the initial data.


Introduction
We consider the evolution of two incompressible fluids with the same viscosity and different densities, moving in a porous two-dimensional medium with constant permeability under the action of gravity according to Darcy's law. After non-dimensionalizing, the equations describing the evolution of density ρ and velocity u are given by (see [11,16] and references therein) We refer to [1,4] for results of the general case. In this note, we assume that at the initial time the two fluids, with densities ρ + and ρ − , are separated by an interface which can be written as the graph of a function over the horizontal axis. That is, Thus, the interface separating the two fluids at the initial time is given by Γ := {(s, z 0 (s))|s ∈ R}. We distinguish the following cases: If ρ + > ρ − , which means that the heavier fluid is on top, we speak of the unstable regime. The case ρ + < ρ − is called the stable regime.
The Muskat problem. Since for given ρ(x, t) at a fixed time t, u is the solution of an elliptic problem by the Biot-Savart law, the Eqs. (1)-(3) describe the evolution of the density in time. Assuming that ρ(x, t) remains in the form (5) for positive times, the system reduces to a non-local evolution problem for the interface Γ , known as the Muskat problem (see [20]). If the sheet can be presented as a graph as above, one can show (see for example [11]) that the equation for z(s, t) is given by The behavior of solutions of (6) depends strongly on the sign of ρ − − ρ + . In the stable case, this equation is locally well-posed in H 3 (R), we refer to [2,11] or [7,9] for an improved regularity. However, in the unstable regime, we have an ill-posed problem, see [6,8,11,21,24]. In particular, in the unstable case, there are no general existence results for (6) known. Thus, the description of (1)-(4) as a free boundary problem seems not suitable for the unstable regime. In [16][17][18] F. Otto used a Lagrangian relaxational approach in the spirit of optimal transportation and he proved the existence of a unique relaxation limitρ in the case of the flat initial datum z 0 ≡ 0. Whilst Otto's relaxation limit does not satisfy the original system (1)-(4) any more, the relaxed density can be thought of as a macroscopic average of an infinitely fine mixture of the two phases ρ ± . More precisely, there is a growing mixing zone around the initial unstable sheet Γ , where the two densities are mixed with average densityρ(x, t) which satisfies an evolution equation (a variant of the 1D Burgers equation). Such a behaviour is reminiscent of the physically expected behaviour in the unstable regime [16,29].
Weak solutions. Taking the curl of (3), we can eliminate the pressure and obtain curl u = −∂ x 1 ρ. This motivates the definition of weak solutions in the following form. Definition 1. Let ρ 0 ∈ L ∞ (R 2 ) and T > 0. We call (ρ, u) ∈ L ∞ (R 2 × [0, T )) a weak solution of (1)- (4) with initial data ρ 0 if In [10] infinitely many weak solutions of (1)-(4) were constructed to any initial datum ρ 0 . The construction is a variant of convex integration, as used in [12,14,25], and in particular the result in [10] is in a sense the IPM-analogue of the Scheffer-Shnirelman construction for the Euler equations [22,23,28]. However, these weak solutions do not retain the geometric structure of initial data of the type (5), and in particular the density ρ may exceed the initial densities ρ ± .
Admissibility and mixing solutions. Motivated by the analogous development of admissible weak solutions for the Euler equations [13] as well as the result of Otto in [16], in [27] admissible weak solutions were introduced by the second author as weak solutions (ρ, u) such that ρ(x, t) ∈ [ρ − , ρ + ] for a.e. (x, t), In [27] the second author showed that there exist infinitely many admissible weak solutions for the Muskat problem in the unstable regime with flat initial data by the convex integration method. Moreover, an interesting connection between such admissible weak solutions and the relaxation limit of Otto is provided by the concept of subsolution (see below in Definition 2). In a nutshell, weak solutions constructed by convex integration arise by adding high-frequency spatially localized perturbations to an initialρ(x, t), which can be thought of as an averaged density: by increasing the frequency of the perturbations, one can easily construct a sequence of admissible weak solutions (ρ k , u k ), such that ρ k * ρ in L ∞ , see [27]. Whereas the construction of weak solutions from strict subsolutions is by now very well understood in the general setting (see for example [25]), constructing strict subsolutions to the initial value problem still requires ad-hoc methods [3,15,26,27] and no general technique seems to exist.
Recently, the result from [27] was generalized to arbitrary initial curves z 0 by Castro, Cordoba, and Faraco in [5]. The main theorem in [5] states that for each z 0 ∈ H 5 (R) and for ρ + > ρ − , there exist infinitely many admissible weak solutions to (1)-(4) with initial data (5). The key point in the proof is to show the existence of certain strict subsolutions, which are in a sense the geometrically nonlinear analogues of the subsolution constructed in [27]: the two ingredients defining the subsolution are a density and an evolving sheet (as in the original Muskat problem), whose translates are then level-sets of the density (see also Sect. 2 below). The density is chosen exactly as in [27], but the evolving sheet needs to solve a nonlinear and nonlocal evolution equation ∂ t z = F(u) (see (1.11)-(1.12) in [5]) -the analysis of this equation is the central part of the proof in [5].
The main result. The aim of this paper is to give an alternative and considerably simpler proof of the main result from [5] for the unstable case. The key difference is that we allow the density to be piecewise constant -in turn, rather than having to prove local well-posedness for a non-linear and non-local evolution equation, it suffices to obtain expressions for the velocity and acceleration of a double-sheet at time t = 0. Our construction is similar in spirit to fan-subsolutions, introduced for flat shock-waves for the compressible Euler equations in [15]. The advantage is not only the considerably shorter proof, but also the lower regularity requirement on the initial curve: we require C 3,α with decay at infinity, in contrast to H 5 in [5]. Finally, we extend our result also to the stable case, provided the initial interface is not horizontal flat (thus, extending the observation made in [5] Section 7 concerning flat non-horizontal interfaces in the stable case).
Our assumption on the initial datum is that the initial interface is asymptotically flat with some given slope β ∈ R, i.e. ρ 0 is given by (5) with for some z 0 with sufficiently fast decay at infinity. More precisely, for any α > 0 define the seminorm and for any k ∈ N the norm We denote by C k,α (i) In the unstable case ρ + > ρ − , for each β ∈ R, there exists T * > 0 such that there exist infinitely many admissible weak solutions to (1)-(4) in [0, T * ). (ii) In the stable case ρ + < ρ − , whenever β = 0 and ∂ sz0 L ∞ < |β| there exists T * > 0 such that there exist infinitely many admissible weak solutions to (1)-(4) in [0, T * ).
As pointed out above, the advantage of our method is the simplicity of the proof and the lower regularity requirement for the initial curve. However, this comes at a small price: as will be explained below in Sect. 2, the admissible weak solutions obtained in our Theorem (as well as those obtained in [5,27]) have the common feature, that there is an expanding mixing zone Ω mi x (t) concentrating on the initial interface Γ 0 at time t = 0, where the two fluids are "infinitely mixed". The rate of expansion of the mixing zone in the unstable case is given by for some c > 0. The constructions in [5,27] admit any c < 2, and indeed, c max = 2 seems the maximal expansion rate possible (see [27] and the discussion following Theorem 2 below). In contrast, our construction admits only c < 1. However, this is only a problem for the simplest possible choice of piecewise constant density in (14) -we will show in Sect. 5 that with a more general piecewise constant density any expansion rate c < 2 is reachable.
We also point out that, just like in [5], our result is local in time: there is a short time of existence [0, T * ]; moreover T * → 0 as we reach the maximal speed c → c max . This is at variance with the result in [27], which is global in time.
The paper is organized as follows. In Sect. 2 we recall the notion of a subsolution for (1)-(4) and show in Theorem 3 that under appropriate estimates for the interface z(s, t) as time t → 0, we are able to construct a subsolution and thus prove our main result Theorem 1. In Sect. 3 we recall the expression for the normal component of the velocity obtained by the Biot-Savart law for piecewise continuous densities, and provide Schauder-type estimates for the associated integral operators. Section 4 is devoted to the construction of the interface curve z(s, t) by evaluating the velocity and symmetrized acceleration at time t = 0. Finally, we generalize these results to more general piecewise constant densities in Sect. 5.

Subsolutions for the IPM Equations
We start by recalling the general strategy for the construction of weak solutions, as followed in [5,27]. The basic idea is to construct a suitable admissible subsolution, and then apply the general machinery of convex integration.
Observe that if (ρ, u) is a solution of (1)-(4), then so is (ρ,ũ) given bỹ Then, by choosing we may assume that the Muskat-type initial datum (5) is given by ρ ± = ±1 (in the stable case the signs are obviously swapped). Under this normalization, admissibility amounts to the requirement |ρ| 1 for a.e. (x, t).
Admissible subsolutions lead to the existence of infinitely many admissible weak solutions by the Baire category method for convex integration, see for instance the Appendix in [27]. As pointed out in [5], a slight modification of the general technique leads to the following statement: Theorem 2. Suppose there exists an admissible subsolution (ρ,ū,m) to (1)- (4). Then there exist infinitely many admissible weak solutions (ρ, u) with the following additional mixing property: For any r > 0, 0 < t 0 < T and x 0 ∈ R 2 such that B : Furthermore, there exists a sequence of such admissible weak solutions (ρ k , u k ) such that ρ k * ρ as k → ∞.
Thus, the crux of the matter is the construction of an admissible subsolution. In [27] the x 1 -invariance of the initial curve z 0 ≡ 0 simplifies the construction of a subsolution. Indeed, assuming that (ρ, u, m) is a function of (x 2 , t) only, the equation ∂ t ρ +div m = 0 together with maximizing the constraint (11) leads to Burger's equation ∂ t ρ +c∂ x 2 1 2 ρ 2 = 0 with 0 c < 2. This equation admits a continuous rarefaction wave solution for the density (12) in the unstable case, whereas in the stable case we merely obtain the stationary shockwave Therefore c can be thought of as a weak notion of mixing speed, in the following sense: the general structure of weak solutions corresponding to this subsolution will be that there are three time-dependent regions: Ω + (t), Ω − (t) and Ω mi x (t), given by The three open sets in Definition 2 are then In Ω ± the density is given by the constant value ρ ± = ±1, and in the mixing zone Ω mi x the two fluids are completely mixed -see Section 4 in [27] and Sections 2-3 in [5]. In the above x 1 -invariant setting from [27] the curve z is simply stationary, i.e. z(s, t) = z 0 (s) ≡ 0. Furthermore, it was shown in [27] that for x 1 -invariant subsolutions c = 2 is the maximal possible speed, and it was conjectured, based on similarities with the Lagrangian relaxation framework of Otto in [16,17] that the maximal mixing speed could be used as a selection criterion.
In [5] the construction from [27] was generalized to non-flat initial curves z 0 whilst retaining the structure (12) in the mixing zone Ω mi x . More precisely, the density is chosen as a linear interpolation between ρ + = 1 and ρ − = −1. In this case, however, z = z(s, t) has to solve a rather complicated evolution equation in time, which arises as a spatial average of the original Muskat evolution kernel. We wish to emphasize that the evolution equation obtained in [5] is not necessarily a canonical choice, but rather arises from the specific ansatz used for ρ -indeed, we show below that simpler choices for the profile ρ reduce the existence of a subsolution to differential inequalities which can be solved by prescribing velocity and acceleration of the interfaces at time t = 0. Indeed, given z : R × [0, T ] → R define Ω ± (t) and Ω mi x (t) as in (13) and set where ρ + = 1, ρ − = −1 in the unstable case, and ρ + = −1, ρ − = 1 in the stable case. This definition of ρ already determines the velocity u by kinematic part of (10), namely the Biot-Savart rule (see Sect. 3 below) Note that ρ is piecewise constant, with jump discontinuities across two interfaces It is well known [11] that, provided the interfaces are sufficiently regular, the solution u to (15) is globally bounded, smooth in R 2 \ (Γ + ∪ Γ − ) with well-defined traces on Γ ± , and the normal component is continuous across the interfaces. In particular, it follows that the normal velocity component is a globally defined bounded and continuous function. In particular we set Our main result in this section is as follows: In the stable case assume in addition that z(·, 0) L ∞ < |β|. Then, there exists T * ∈ (0, T ] such that there exists an admissible subsolution for (1)-(4) on [0, T * ) with initial datum ρ 0 given by (5) with z 0 = z| t=0 . Furthermore, the density of the subsolution can be chosen to satisfy (13)- (14) with any 0 < c < c max , where Remark 1. We note that the time of existence T * > 0 depends on c and in particular Proof. Given z = z(s, t), c > 0 and ρ(x, t) (defined by (13)- (14)), the velocity u is determined by (15). Therefore it remains to define m so that (10)- (11) are satisfied in for some γ = γ (x, t), with γ ≡ 0 in Ω ± . Then (11) amounts to the condition whereas (10) is equivalent to div γ = 0 in Ω mi x together with two jump conditions where [·] Γ ± denotes the jump on Γ ± . Noting that u ν in (17) is globally well-defined and continuous, the jump conditions become where is the tangential derivative of g along curves defined by z. We treat the unstable and stable cases separately.
and observe that In order to satisfy (22) we then set and, more generally, for λ ∈ [−ct, ct] as t → 0. Therefore, for any 0 < c < 1 we deduce that This concludes the proof in the unstable case. Stable case. Define the one-parameter family of diffeomorphisms and observe that In order to satisfy (23) we set and, more generally, for λ ∈ [−ct, ct] Then From the assumptions (18)- (19) it follows that Therefore we obtain we have |∇g| < 1/2 for sufficiently small t > 0. This concludes the proof in the stable case.

The Velocity u
In this section we derive a concrete representation formula for the velocity u and for the normal velocity component u ν defined in (17), where u is the solution of the system (15). It is well-known [19] that for sufficiently smooth ρ the solution v of can be written using the Biot-Savart kernel as If the density ρ is piecewise constant, with a jump across a sufficiently smooth interface Γ , the expression for v(x) for x / ∈ Γ can be derived by formally writing ∂ x 1 ρ as a delta distribution supported on Γ [11]. More precisely, in [11] (see Section 2 therein) the following expression was derived for the normal velocity component v ν under the assumption that the interface is given by a graph Γ = {(s, z(s)) : s ∈ R} with z ∈ C 1,α (R) and For any where the principal value refers to the limit lim R→∞ R −R . For the convenience of the reader we recall the argument leading up to formula (26). First of all, by writing ∂ x 1 ρ as a delta distribution supported on Γ , from (25) for all x / ∈ Γ . Then, by using that Then, for the density ρ = ρ(x, t) defined in (14) with interfaces Γ ± in (16), the normal velocity component (17) on Γ ± is given by the expression where Motivated by this expression, consider the following Φ-weighted variant of the Hilbert transform: T for some bounded weight-function Φ = Φ(ξ, s). For the weight we use the following norms: first of all we assume that Φ ∞ (s) We introduce the norms where we use the convention that Φ(ξ, ·) denotes a norm in the second argument only and Φ denotes a norm joint in both variables. In particular the Hölder-continuity of ∂ k s Φ in both variables ξ, s is required in the norm |||Φ||| k,α . Accordingly, we define the spaces We have the following version of the classical estimate on the Hilbert transform T 1 = H∇ on Hölder-spaces: Moreover, for any k ∈ N, f ∈ C k+1,α * where the constant depends only on k and α.
Proof. We start the proof by rewriting the principal value integral in (28) as a sum of absolutely convergent terms. To this end we split the integral R dξ = |ξ |<1 dξ + |ξ |>1 dξ and integrate by parts in the second term. We obtain whereΦ andΦ are related to Φ as in (29). It is easy to see that For the term I 4 observe that if |ξ | < 1 2 |s| or |ξ | > 3 2 |s|, then |s−ξ | > 1 2 |s|. Therefore we have for any |s| > 2 This concludes the proof of (30). The Hölder continuity of I 2 , I 3 and I 4 is easily handled analogously and leads to the estimates Next, we consider I 1 = I 1 (s) and write for simplicity g(s) = ∂ s f (s). For |η| < 1/2 let s = s − η and write where I 11 is an integral over intervals of total length ∼ |η| on which |s − ξ | > 1/2. Therefore Next, we write, with r = 2|η| We can estimate each term as follows: We conclude that and this finally proves (31) for k = 0. For k 1 the estimate follows from differentiating the terms I 1 (s), . . . , I 4 (s) with respect to s and applying the Leibniz rule.

Construction of the Curve z
In this section we construct a function z = z(s, t) satisfying the conditions of Theorem 3. In order to motivate the construction, observe that (18)- (19) suggest that it suffices to specify z up to order t 2 . Therefore we start by formally calculating the expressions for the initial velocity u ± ν | t=0 and initial symmetrized acceleration ∂ t | t=0 for some β ∈ R and α ∈ (0, 1). Using the expression (27) and the notation introduced in (28) we have where ρ ± = ±1 in the unstable case and ρ ± = ∓1 in the stable case. This difference in sign has no effect on the computations and on Theorem 5 below, therefore we will from now on treat the unstable case without loss of generality. Hence, in particular and z 0 (s) = z(s, 0). Observe that although Φ ± (ξ, s) → Φ 0 (ξ, s) as t → 0 for any ξ = 0, the limit is not uniform in ξ , therefore in particular Φ ± Φ 0 in the norm of W 0 . Nevertheless we have Lemma 2. Assume that z(s, t) = βs +z(s, t) withz ∈ C 0 ([0, T ); C 1,α * (R)) for some β ∈ R and α ∈ (0, 1). Then for any f ∈ C 1,α * (R) Proof. In analogy with (32) we set . In the following we consider without loss of generality Φ + (t) − Φ 0 .
As pointed out above, the limit lim t→0 Φ + (t) is not uniform in ξ because of the singularity at ξ = 0. Therefore we need to modify the argument in the proof of (30). To this end recall the decomposition in the proof of Theorem 4 and focus for the moment on the term Using the definition of Φ + and Φ 0 we write for ξ = 0 It is easy to see that sup ξ,s,t |Φ + | 2. Moreover sup ξ,s On the other hand, by splitting the integral |ξ |<1 dξ = |ξ |<t 1/2 dξ + t 1/2 <|ξ |<1 dξ we can estimate Next, from (35) we equally deduce Concerning I 2 , note that for |ξ | > 1 Since also sup |ξ |>1,s∈R Finally, let us look at I 4 , which requires bounding ξ∂ ξ (Φ + − Φ 0 ). Observe that By a simple calculation it follows lim t→0 sup |ξ |>1,s∈R |ξ∂ ξ (Φ + − Φ 0 )(ξ, s, t)| = 0 and hence This concludes the proof of the Lemma.
Next, in order to evaluate ∂ ∂t t=0 where z 0 (s) = ∂ t z(s, 0) and Z 0 (ξ, s) = z 0 (s)−z 0 (s−ξ) ξ . For the regular part Δ reg Φ we have and α ∈ (0, 1). Then for any f ∈ C 1,α and by assumption is uniformly bounded in s, ξ as t → 0, the first summand can be dealt with exactly as in the proof of Lemma 2. On the other hand the second summand converges to zero uniformly in s, ξ as t → 0. This concludes the proof.
Next we analyse the singular part Δ sing Φ: for some β ∈ R and α ∈ (0, 1). Then, for any f ∈ C 2 (R) for any s ∈ R, where Proof. We begin by performing the change of variables ξ → ξ t in the integral: .
Since sup ξ,s,t |Z t | ∂ s z L ∞ and where the constant depends only on ∂ s z L ∞ and on c > 0. Furthermore, since we deduce that for any ξ, s ∈ R so that, from the Lebesgue dominated convergence theorem we deduce that Noting that the bound (39) applies also to Ψ 0 and ∂ s z 0 is independent of ξ , we may evaluate the integral R Ψ 0 dξ using elementary methods. Indeed, we calculate for any constants a ∈ R and c > 0 hence the proof of (37) follows.
In order to show (38) we again start with the decomposition 1 t T Δ sing Φ f = I 1 + I 2 + I 3 + I 4 as in the proof of Theorem 4. We claim that and for k = 2, 3, 4 sup The proof of (41) follows from the observation that 1 t Δ sing Φ and ξ∂ ξ 1 t Δ sing Φ are bounded uniformly in s ∈ R, |ξ | 1. The claim (40) is equivalent to showing lim t→0 J t = 0, where Let ε > 0 and fix R > 1 so that |ξ |>R |Ψ t | dξ < ε for all t 0 (by the bound (39) this is possible). Moreover, fix 0 < δ < 1 so that for all 0 < |η| < δ. This is possible if f ∈ C 2,α * (R) since we can write Then, for any t > 0 with t R < δ we have Using once more the bound (39) we deduce that lim sup t→0 J t C(1 + f * 2,α )ε. Since ε > 0 was arbitrary, it follows that lim t→0 J t = 0, concluding the proof of (40) and thence the proof of (38).
Next, we calculate: Note that, since ∂ s z 0 ∈ C 2,α (R), the function Z 0 (ξ, s) satisfies Z 0 ∈ C 2,α (R 2 ) and this implies Φ 0 ∈ C 2,α (R 2 ). Furthermore, proceeding as above, it follows easily that ξ∂ k s Z 0 is bounded uniformly for s ∈ R, |ξ | > 1 for k = 1, 2, and hence the same is true for ξ∂ k s Φ 0 . Analogously, ξ∂ ξ ∂ k s Z 0 is bounded uniformly for s ∈ R, |ξ | > 1 for k = 0, 1, 2, hence the same is true for ξ∂ ξ ∂ k s Φ 0 . This implies that ∂ s Φ 0 , ∂ 2 s Φ 0 ∈ W 0 . In the same manner we deduce that for k = 1, 2 |ξ ||∂ k We next proceed with Φ 1 . First of all, from Theorem 4 and the above we deduce that z 1 ∈ C 2,α * . Comparing the expressions for Φ 1 and ∂ s Φ 0 we see that the estimates for Φ 1 may be obtained exactly as above, by replacing ∂ s z 0 by z 1 where appropriate and noting that both functions are in C 2,α . In this way we deduce that Φ 1 ∈ W 1,α .

Symmetric Piecewise Constant Densities
The subsolution (and corresponding admissible weak solutions) constructed in the previous sections have a mixing zone with a maximal expansion rate of c max = 1 in the unstable case, whereas the maximal expansion rate reachable with the construction in [5] is c max = 2. In this section we show that with a more general piecewise constant density the method of this paper is applicable to reach the expansion rate c max = 2 as well -indeed, this is easily achieved by approximating the linear density function from [5] by a piecewise constant density. From now on we restrict attention to the unstable case, with ρ + = 1 and ρ − = −1.
Let N ∈ N and define 2N interfaces where 0 < c 1 < c 2 < . . . < c N are arbitrary velocities and z(s, t) is the parametrization of a curve for t ∈ [0, T ]. The open regions between neighbouring interfaces are defined as and Analogously to (14) we set so that we have a constant density jump of 1 N across each boundary Γ ±i . The mixing zone is then With the density defined in this way, the velocity u can be obtained as in Sect. 3, in particular we have the expression for the normal velocity component: 2 and we have set c −i = −c i for i = 1, . . . , N . Theorem 6. Suppose that z(s, t) = βs +z(s, t) withz ∈ C 1 ([0, T ); C 1,α * (R)) satisfies Moreover, let Then there exists T * ∈ (0, T ) such that there exists an admissible subsolution for (1)-(4) on [0, T * ] with initial datum ρ 0 given by (5) in the unstable case with z 0 = z| t=0 . Furthermore, the density of the subsolution satisfies (43)-(45).
This concludes the proof.
Finally, observe that for the subsolution obtained in Theorem 6 the rate of expansion of the mixing zone is given by c n < 2n−1 n = 2 − 1 n , so that any expansion rate c < 2 is obtainable by choosing n sufficiently large.