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 $C^{3,\alpha}$-regularity of the interface in the unstable regime and for all non-horizontal data with $C^{3,\alpha}$-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 [10,5] and references therein) u`∇p "´p0, ρq , ρpx, 0q " ρ 0 pxq . (4) 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 Γ :" tps, z 0 psqq|s P Ru. 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 ρpx, tq at a fixed time t, u is the solution of an elliptic problem by the Biot-Savart law, the equations (1)- (3) describe the evolution of the density in time. Assuming that ρpx, tq 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. If the sheet can Date: September 18, 2017. 1 be presented as a graph as above, one can show (see for example [5]) that the equation for zps, tq is given by B t zps, tq " ρ´´ρ2 π ż 8 8 pB s zps, tq´B s zpξ, tqqps´ξq ps´ξq 2`p zps, tq´zpξ, tqq 2 dξ.
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 pRq, see [5] or [3] for an improved regularity. However, in the unstable regime, we have an ill-posed problem, see [13,5]. 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 [10,11] 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ρpx, tq 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 [19,10].
In [4] 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 [6,8], and in particular the result in [4] is in a sense the IPM-analogue of the Scheffer-Shnirelman construction for the Euler equations [14,15,18]. 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 [7] as well as the result of Otto in [10], in [17] admissible weak solutions were introduced by the second author as weak solutions pρ, uq such that (7) ρpx, tq P rρ´, ρ`s for a.e. px, tq, (or ρ P rρ`, ρ´s in the stable case, where ρ`ă ρ´).
In [17] 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.1). In a nutshell, weak solutions constructed by convex integration arise by adding high-frequency spatially localized perturbations to an initial ρpx, tq, 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 pρ k , u k q, such that ρ ká ρ in L 8 , see [17]. Whereas the construction of weak solutions from strict subsolutions is by now very well understood in the general setting, constructing strict subsolutions to the initial value problem still requires ad-hoc methods [16,17,1,9] and no general technique seems to exist.
Recently, the result from [17] was generalized to arbitrary initial curves z 0 by Castro, Cordoba and Faraco in [2]. The main theorem in [2] states that for each z 0 P H 5 pRq 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 [17]: 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 Section 2 below). The density is chosen exactly as in [17], but the evolving sheet needs to solve a nonlinear and nonlocal evolution equation B t z " Fpuq (see (1.11)-(1.12) in [2]) -the analysis of this equation is the central part of the proof in [2].
The main result. The aim of this paper is to give an alternative and considerably simpler proof of the main result from [2] for the unstable case.
The key difference is that we allow the density to be piecewise constantin 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 [9]. 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 [9]) . Finally, we extend our result also to the stable case, provided the initial interface is not horizontal flat (thus, extending the observation made in [2] Section 7 concerning flat non-horizontal interfaces in the stable case).
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 Section 2, the admissible weak solutions obtained in our Theorem (as well as those obtained in [17] and [2]) have the common feature, that there is an expanding mixing zone Ω mix ptq 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 ρ`´ρ2 c for some c ą 0. The constructions in [17] and [2] admit any c ă 2, and indeed, c max " 2 seems the maximal expansion rate possible (see [17] and the discussion following Theorem 2.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 Section 5 that with a more general piecewise constant density any expansion rate c ă 2 is reachable.
We also point out that, just like in [2], our result is local in time: there is a short time of existence r0, T˚s; moreover T˚Ñ 0 as we reach the maximal speed c Ñ c max . This is at variance with the result in [17], which is global in time.
The paper is organized as follows. In Section 2 we recall the notion of a subsolution for (1)-(4) and show in Theorem 2.3 that under appropriate estimates for the interface zps, tq as time t Ñ 0, we are able to construct a subsolution and thus prove our main result Theorem 1.2. In Section 3 we recall the expression for the normal component of the velocity obtained by the Biot-Savart law for piecewise continuous densities, and provide Schaudertype estimates for the associated integral operators. Section 4 is devoted to the construction of the interface curve zps, tq by evaluating the velocity and symmetrized acceleration at time t " 0. Finally we generalize these results to more general piecewise constant densities in Section 5.

Subsolutions for the IPM Equations
We start by recalling the general strategy for the construction of weak solutions, as followed in [2] and [17]. The basic idea is to construct a suitable admissible subsolution, and then apply the general machinery of convex integration.
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. px, tq.
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 [17]. As pointed out in [2], a slight modification of the general technique leads to the following statement: Theorem 2.2. Suppose there exists an admissible subsolution pρ,ū,mq to (1)- (4) . Then there exist infinitely many admissible weak solutions pρ, uq with the following additional mixing property: For any r ą 0, 0 ă t 0 ă T and x 0 P R 2 such that B :" B r px 0 , t 0 q Ă Ω mix , both sets tpx, tq P B : ρpx, tq "˘1u have strictly positive Lebesgue measure.
Furthermore, there exists a sequence of such admissible weak solutions pρ k , u k q such that ρ káρ as k Ñ 8.
Thus, the crux of the matter is the construction of an admissible subsolution. In [17] the x 1 -invariance of the initial curve z 0 " 0 simplifies the construction of a subsolution. Indeed, assuming that pρ, u, mq is a function of px 2 , tq only, the equation B t ρ`div m " 0 together with maximizing the constraint (11) leads to Burger's equation B t ρ`cB x 2 1 2 ρ 2 " 0 with 0 ď c ă 2. This equation admits a continuous rarefaction wave solution for the density (12) ρpx, tq " in the unstable case, whereas in the stable case we merely obtain the stationary shock-wave 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: Ω`ptq, Ω´ptq and Ω mix ptq, given by for some curve zp¨, tq, with Ω mix ptq expanding with speed c. The three open sets in Definition 2.1 are then In Ω˘the density is given by the constant value ρ˘"˘1, and in the mixing zone Ω mix the two fluids are completely mixed -see Section 4 in [17] and Sections 2-3 in [2]. In the above x 1 -invariant setting from [17] the curve z is simply stationary, i.e. zps, tq " z 0 psq " 0. Furthermore it was shown in [17] 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 [10,11] that the maximal mixing speed could be used as a selection criterion.
In [2] the construction from [17] was generalized to non-flat initial curves z 0 whilst retaining the structure (12) in the mixing zone Ω mix . More precisely, the density is chosen as a linear interpolation between ρ`" 1 and ρ´"´1. In this case, however, z " zps, tq 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 [2] 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ˆr0, T s Ñ R define Ω˘ptq and Ω mix ptq as in (13) and set ρpx, tq " 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 Section 3 below) div u " 0, curl u "´B x 1 ρ. (15) Note that ρ is piecewise constant, with jump discontinuities across two interfaces (16) Γ˘ptq " tps, zps, tq˘ctq : s P Ru .
It is well known [5] that, provided the interfaces are sufficiently regular, the solution u to (15) is globally bounded, smooth in R 2 zpΓ`Y Γ´q with welldefined traces on Γ˘, and the normal component is continuous across the interfaces. In particular it follows that the normal velocity component s a globally defined bounded and continuous function. In particular we set uν " u ν | Γ˘, i.e. uν ps, tq :" u ν ps, zps, tq˘ct, tq. Our main result in this section is as follows: In the stable case assume in addition that }zp¨, 0q} L 8 ă |β|. Then, there exists T˚P p0, T s such that there exists an admissible subsolution for (1)-(4) on r0, T˚q 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 in the stable case.
Remark 2.4. We note that the time of existence T˚ą 0 depends on c and in particular T˚Ñ 0 as c Ñ c max .
Proof. Given z " zps, tq, c ą 0 and ρpx, tq (defined by (13)- (14)), the velocity u is determined by (15). Therefore it remains to define m so that (10)- (11) are satisfied in p0, T˚qˆR 2 , with (11) a strict inequality in Ω mix . Set m " ρu´p1´ρ 2 qpγ`1 2 e 2 q for some γ " γpx, tq, with γ " 0 in Ω˘. Then (11) amounts to the condition whereas (10) is equivalent to div γ " 0 in Ω mix together with two jump conditions where r¨s Γ˘d enotes the jump on Γ˘. Noting that u ν in (17) is globally well-defined and continuous, the jump conditions become and γν denotes the one-sided limit lim Ω mix Qx 1 Ñx γ ν pxq for x P Γ˘. Choosing g˙f or some function g P C 1 pΩ mix q and noting that is the tangential derivative of g along curves defined by z. We treat the unstable and stable cases separately.
This concludes the proof in the unstable case. Stable case. Define the one-parameter family of diffeomorphisms Moreover, since Φ t ps, λq P tx 2 " zpx 1 q`λu for any s P R, it follows that Φ t maps Rˆr´ct, cts onto Ω mix . Set gps, λ, tq " gpΦ t pxq, tq.

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 [12] that for sufficiently smooth ρ the solution v of can be written using the Biot-Savart kernel as vpxq :" BSp´B x 1 ρq :" 1 2π If the density ρ is piecewise constant, with a jump across a sufficiently smooth interface Γ, the expression for vpxq for x R Γ can be derived by formally writing B x 1 ρ as a delta distribution supported on Γ [5]. More precisely, in [5](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 Γ " tps, zpsqq : s P Ru with z P C 1,α pRq and For any x " px 1 , x 2 q P R 2 we have v ν pxq :" vpxq¨ˆ´B 1 zpx 1 q 1" where the principal value refers to the limit lim RÑ8 ş Ŕ R .
For the convenience of the reader we recall the argument leading up to formula (26). First of all, by writing B x 1 ρ as a delta distribution supported on Γ, from (25) one obtains for all x R Γ. Then, by using that from which (26) follows. Then, for the density ρ " ρpx, tq defined in (14) with interfaces Γ˘in (16), the normal velocity component (17)  where Φ˘pξ, s, tq " ξ 2 ξ 2`p zps´ξ, tq´zps, tqq 2`ξ 2 ξ 2`p zps´ξ, tq´zps, tq¯2ctq 2 . Motivated by this expression, consider the following Φ-weighted variant of the Hilbert transform: for some bounded weight-function Φ " Φpξ, sq. For the weight we use the following norms: first of all we assume that Φ 8 psq :" lim |ξ|Ñ8 Φpξ, sq exists, Φp¨, sq P C 1 pRzt0uq, and set We introduce the norms where we use the convention that }Φpξ,¨q} denotes a norm in the second argument only and }Φ} denotes a norm joint in both variables. In particular the Hölder-continuity of B k s Φ in both variables ξ, s is required in the norm Φ~k ,α . Accordingly, we define the spaces W 0 " tΦ P L 8 pR 2 q : Φ 8 and B ξ Φ exist, with~Φ~0 ă 8u, W k,α " tΦ P W 0 :~Φ~k ,α ă 8u We have the following version of the classical estimate on the Hilbert transform T 1 " H∇ on Hölder-spaces: Theorem 3.1. For any α ą 0, f P C 1,α pRq and Φ P W 0 we have Moreover, for any k P N, f P C k`1,α pRq and Φ P W k,α 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).
This concludes the proof of (30).
Next, we write, with r " 2|η| We can estimate each term as follows: We conclude that }I 1 }α ď C}Φ} α }B s f }α, and this finally proves (31) for k " 0. For k ě 1 the estimate follows from differentiating the terms I 1 psq, . . . , I 4 psq with respect to s and applying the Leibniz rule.

Construction of the curve z
In this section we construct a function z " zps, tq satisfying the conditions of Theorem 2.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 B t | t"0 uν`uν 2 . Let z " zps, tq " βs`zps, tq withz P C 2 pr0, T q; C 1,α pRqq for some β P R and α P p0, 1q. 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 4.4 below, therefore we will from now on treat the unstable case without loss of generality. Hence, in particular Φ 0 pξ, sq :" 2ξ 2 ξ 2`p z 0 ps´ξq´z 0 psqq 2 and z 0 psq " zps, 0q. Observe that although Φ˘pξ, sq Ñ Φ 0 pξ, sq 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 4.1. Assume that zps, tq " βs`zps, tq withz P C 0 pr0, T q; C 1,α pRqq for some β P R and α P p0, 1q. Then for any f P C 1,α pRq Proof. In analogy with (32) we set . In the following we consider without loss of generality Φ`ptq´Φ 0 .
As pointed out above, the limit lim tÑ0 Φ`ptq 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 3.1 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 On the other hand, by splitting the integral ş |ξ|ă1 dξ " Hence lim tÑ0 sup s p1`|s| 1`α q|I 1 ps, tq| " 0.

Now write
0 q is uniformly bounded in s, ξ as t Ñ 0, the first summand can be dealt with exactly as in the proof of Lemma 4.1. On the other hand the second summand converges to zero uniformly in s, ξ as t Ñ 0. This concludes the proof.

Proof.
Step 1: Estimating z 1 and z 2 . We begin by showing that z 1 P C 2,α pRq and z 2 P C 1,α pRq. First of all, since B s z 0 P C 2 pRq with B 2 s z 0 P C 1,α pRq, it follows that σ P C 2 pRq with σ, B s σ, B 2 s σ P L 8 pRq and consequently cσB 2 s z 0 P C 1,α pRq. Therefore, according to estimate (31) in Theorem 3.1 it suffices to show that Φ 0 P W 2,α and Φ 1 P W 1,α , where Φ 0 and Φ 1 are defined in (33) and (36) above.
is bounded uniformly in ξ, s. Similarly, we calculate We deduce that alsoΦ 0 is uniformly bounded. This shows that Φ 0 P W 0 .
Next, we calculate: Note that, since B s z 0 P C 2,α pRq, the function Z 0 pξ, sq satisfies Z 0 P C 2,α pR 2 q and this implies Φ 0 P C 2,α pR 2 q. Furthermore, proceeding as above, it follows easily that ξB k s Z 0 is bounded uniformly for s P R, |ξ| ą 1 for k " 1, 2, and hence the same is true for ξB k s Φ 0 . Analogously, ξB ξ B k s Z 0 is bounded uniformly for s P R, |ξ| ą 1 for k " 0, 1, 2, hence the same is true for In the same manner we deduce that for k " 1, 2 |ξ||B k s Z 0 pξ, s 1 q´B k s Z 0 pξ, s 2 q| ď|B k s z 0 ps 1 q´B k s z 0 ps 2 q| |B k s z 0 ps 1´ξ q´B k s z 0 ps 2´ξ q| ď2rB k s z 0 s α |s 1´s2 | α and similarly, for |ξ| ą 1 and k " 1, 2 We next proceed with Φ 1 . First of all, from Theorem 3.1 and the above we deduce that z 1 P C 2,α . Comparing the expressions for Φ 1 and B s Φ 0 we see that the estimates for Φ 1 may be obtained exactly as above, by replacing B s z 0 by z 1 where appropriate and noting that both functions are in C 2,α . In this way we deduce that Φ 1 P 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 [2] 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 [2] by a piecewise constant density. From now on we restrict attention to the unstable case, with ρ`" 1 and ρ´"´1.
Analogously to (14) we set (45) ρpx, tq " i N for x P Ω i ptq, i "´N, . . . , N, 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 Section 3, in particular we have the expression for the normal velocity component: Φ i,j pξ, s, tq " ξ 2 ξ 2`p zps´ξ, tq´zps, tq´pc i´cj qtq 2 ξ 2 ξ 2`p zps´ξ, tq´zps, tq´pc i`cj qtq 2 and we have set c´i "´c i for i " 1, . . . , N .
Proof. We proceed as in the proof of Theorem 2.3 and set m " ρu´p1´ρ 2 qpγ`1 2 e 2 q, where g pN q " g p´N q " 0 and g piq P C 1 pΩ i q for i "´pN´1q . . . pN´1q are to be determined. Then, (11) amounts to the conditions and (10) reduces to jump conditions (20) on each interface: for any i " 1, . . . , N with h p˘iq ps, tq " and B τ gpx, tq " B x 1 gpx, tq`B x 2 gpx, tqB x 1 zpx 1 , tq.
This concludes the proof.
The proof is entirely analogous to the proof of Theorem 4.4, based on the above calculations and Lemmas 4.1, 4.2 and 4.3.
Finally, observe that for the subsolution obtained in Theorem 5.1 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.