Weak Solutions to the Muskat Problem with Surface Tension Via Optimal Transport

Inspired by recent works on the threshold dynamics scheme for multi-phase mean curvature flow (by Esedoḡlu–Otto and Laux–Otto), we introduce a novel framework to approximate solutions of the Muskat problem with surface tension. Our approach is based on interpreting the Muskat problem as a gradient flow in a product Wasserstein space. This perspective allows us to construct weak solutions via a minimizing movements scheme. Rather than working directly with the singular surface tension force, we instead relax the perimeter functional with the heat content energy approximation of Esedoḡlu–Otto. The heat content energy allows us to show the convergence of the associated minimizing movement scheme in the Wasserstein space, and makes the scheme far more tractable for numerical simulations. Under a typical energy convergence assumption, we show that our scheme converges to weak solutions of the Muskat problem with surface tension. We then conclude the paper with a discussion on some numerical experiments and on equilibrium configurations.


Introduction
The Muskat problem was first introduced by Morris Muskat [27] as a model for the flow of two immiscible fluids through a porous medium.Since its introduction, the Muskat problem has received sustained attention in a variety of fields.The Muskat problem is used to model flows in oil reservoirs (water is injected into the oil well to drive oil extraction), and in hydrology to model flows of groundwater through aquifers.
In this paper we are interested in obtaining the global existence of weak solutions for the Muskat problem with surface tension, based on its gradient flow structure.We begin by introducing a variational formulation of the problem, which will motivate our subsequent analysis.The fluid evolution can be written as Darcy's law (1) v i + b −1 i ∇δ ρi E(ρ) = 0, coupled with the continuity equation ( 2) where v i is the velocity of phase i, ρ = (ρ 1 , ρ 2 ) is the collection of relative concentrations for each phase, ∇δ ρi E denotes the gradient of the first variation of the internal energy with respect to ρ i , and b i > 0 (i = 1, 2) denotes constant mobilities.For convenience, throughout the rest of the paper, we will refer to ρ as a collection of density functions; however, one should note that ρ only encodes information about the volume occupied by the fluids and nothing about their mass.
The physical setting for our problem is a bounded, convex open domain Ω ⊂ R d with smooth boundary.We shall suppose that the two fluids fill the entire domain, and that they are confined to Ω for all time.We then take the internal energy to be a sum of three distinct terms: (3) The first term in the energy, describing incompressibility and containment of the fluids, is given by (4) E p (ρ) = 0, if ρ 1 (x) + ρ 2 (x) = 1 for a.e.x ∈ Ω +∞, otherwise.
The immiscibility of the fluids and the surface tension force arise from the highly non-convex interaction energy (5) E s (ρ) = where |Dρ i |(Ω) denotes the total variation of ρ i in Ω and σ > 0 is a surface tension constant.Finally, Φ(ρ) denotes the potential energy of the fluid configuration, i.e.
where Φ 1 , Φ 2 : Ω → R are given Lipschitz continuous potentials.A typical example is when one assumes these to be gravitational potentials, i.e. ( where g i 's are proportional to the specific gravity of each fluid.Although the internal energy is singular, when ρ 1 and ρ 2 are separated by a smooth interface Γ := ∂{ρ 1 > 0} ∩ ∂{ρ 2 > 0}, one can formulate a classical solution to the Muskat problem equations (1-2).In the classical solution, the flow is driven by the pressure variables p i for each phase, which are Lagrange multipliers generated by E p above.The continuity equation becomes 0 and the pressure is determined by solving the free boundary problem where κ denotes the mean curvature of Γ, oriented to be positive when {ρ 2 > 0} is convex at the point, and n denotes the outer normal along ∂Ω and along Γ. Problems like (MP 2 ) received a lot of attention in the past decades.Most of the works focus on the zero surface tension model (σ = 0) and well-posedness of regular solutions with graph property [1,5,7,8,9].In the presence of surface tension, the problem has stronger regularity properties in stable settings [30], but still, topological singularities can occur in finite time, for instance when heavier fluid is placed on top of the lighter one [13].Thus, our aim is to construct global-in-time weak solutions to the Muskat problem (MP 2 ), which exist past the formation of singularities.
To construct global-in-time solutions, we exploit the gradient flow structure of the Muskat problem.As noted by Otto in [28,29], Darcy's law can be approximated by the Euler-Lagrange equation for the minimizing movements scheme (or JKO [17] scheme) with time step size τ > 0, (7) ρ n+1 = arg min where W 2 (ρ i , ρ n i ) denotes the 2-Wasserstein or 2-Monge-Kantorovich distance.In this context, the squared W 2 distance has a physical interpretation as the energy dissipated by friction as the fluids flow through the porous media.
Let us note that Wasserstein gradient flows of energies involving total variation terms have been considered before in the literature, though only in the case of one phase models (see e.g.[24,4]), and hence with no incompressibility or interaction constraints.As a result, the techniques developed in those papers do not appear to be applicable here -the constrained two phase setting adds many additional difficulties.
Indeed, it is not easy to obtain a complete characterization of the solutions to the minimizing movements problem (7).The interaction energy (5) is sufficiently non-convex that problem (7) is non-convex for any τ > 0. As a result, one must be careful in using duality to introduce the pressure as a Lagrange multiplier.Furthermore, we are interested in developing a scheme which could be used for numerical implementations.The formulation (7) is poorly suited for numerical methods.Optimizing over the nonconvex constraint set {ρ 1 , ρ 2 ∈ BV (Ω; {0, 1}) : ρ 1 (x)ρ 2 (x) = 0 a.e.} is extremely difficult.For these reasons, we instead consider a relaxed version of minimizing movements scheme inspired by [12] and [20].

Approximation of the perimeter by the Heat Content
In our analysis, we replace the interaction energy in E s in (7) by the heat content energy Here G ε : R d → R stands for the standard heat kernel (with mean 0 and variance ε > 0) and the densities ρ i are assumed to be defined on all of R d by extending them to zero off of Ω.
The heat content energy was first introduced by Esedoḡlu and Otto in [12].It was observed in [12] that the threshold dynamics, a well known numerical scheme for mean curvature motion introduced by Merriman, Bence, and Osher [26], is precisely a minimizing movements scheme for the heat content energy.Esedoḡlu and Otto also showed that HC ε Γ-converges (with respect to the L 1 topology) to E s as ε → 0. Building off of these results, Laux and Otto showed that under an energy convergence assumption, the threshold dynamics scheme produces weak solutions to the multi-phase motion by mean curvature in the limit ε → 0 [20].
Our goal is to consider such a framework in the context of the Muskat problem by studying the minimizing movements scheme ( 9) where we used the notation As we alluded above, the scheme ( 9) has a number of numerical advantages over (7).Unlike E s which is neither convex nor concave, the heat content is a strictly concave functional of the densities.This concavity can be exploited to simplify numerical implementations, along the same lines as the linearization trick noted in [12].Furthermore, the concavity of the heat content is controlled by ε, thus, if one chooses τ sufficiently small the overall problem becomes displacement convex (c.f.Lemma 2.1).Indeed, displacement convexity is crucial to guarantee that the problem can be solved efficiently with standard numerical techniques.
Although the heat content in principle allows mixing of the phases, we shall show that the discrete in time solutions constructed by the JKO scheme always stay unmixed with a sharp interface between the phases for all time (see Proposition 2.3 below).This phenomenon is due to the fact that HC ε behaves like a strictly concave functional (see Lemma 2.2).Thus, one retains the essential properties of the Muskat problem evolution.
In the context of the Muskat problem, the heat content also has a natural physical interpretation.In a discrete statistical mechanics model with N particles, surface tension can be seen to arise from short range interactions between particles in different phases (cf.[16]).Typically, the discrete surface tension takes the form 1 where P r is the particle index in each phase r, V is some decreasing function, and x i is the location of particle i ( [16]).By taking the limit N → ∞ in above formula, one obtains an analogue of the heat content energy where the kernel G ε is replaced with V .Here it is worth noting that we choose to work with the heat kernel for computational convenience, indeed a different choice of kernel may be more physically relevant.Finally, let us also emphasize that the heat content approach can be naturally extended to the multiphase Muskat Problem evolution (with any number of phases), without incurring any additional difficulties.This includes scenarios where the surface tension force depends on the phases that are interacting [12].Although we do not pursue the multiphase case further in this paper, this presents an interesting direction for future work.

Statement of our main results
From the relaxed minimizing movements scheme (9), we obtain a sequence of discrete in time approximations to the Muskat flow.When we take the time step τ and the heat content approximation parameter ε to zero together, we hope to recover weak solutions to the Muskat problem.Our main results show that this is indeed the case under the assumption that there is no loss of perimeter when passing from discrete to continuous solutions.For the precise statements of the convergence results we refer to Theorem 3.1 and Theorem 3.2.
• Finally, under the assumption of energy convergence Here formally v i corresponds to −b −1 i (∇p + ∇Φ i ).It is challenging to study qualitative or geometric properties of our solutions.We will illustrate heuristically in Section 4 that, even for global minimizers of the energy, there are diverse possibilities depending on the values of the specific gravity and volumes of the two phases.

Remarks on our results
(1) Let us underline the fact that our discrete-time scheme (9) produces minimizers that are characteristic functions of a partition of Ω.To prove this fact, we exploit the strict concavity of the heat content along the admissible set.This seems to be an interesting property in its own right, and ensures that numerical implementations of the scheme maintain a sharp interface at every time step.(2) From the point of view of Wasserstein gradient flows, an interesting remark on our results is that while one cannot expect strong compactness for the interpolated densities (ρ ε 1 , ρ ε 2 ) in the case when ε > 0 is fixed and τ ↓ 0, when sending both parameters to 0 in the same time, we regain the strong compactness.This phenomenon is mainly due to the fact that in the limit as max{τ, ε} ↓ 0 we recover total variation estimates on the densities.This compactness is obtained via a careful Aubin-Lions type argument.
(3) The energy convergence assumption (EC) is rather natural and the same as the ones given in [20] and [21].This assumption ensures that there is no sudden loss of boundary between phases in the limit ε → 0. If there is a loss of interface then one cannot obtain weak solutions.Indeed, if two components of the support of ρ 1 merge into each other and remove a sizable part of its boundary, one can expect a discontinuous change of the pressure term p in the entire domain Ω, creating an inconsistency between the discrete and limiting evolutions.
Replacing the assumption with a direct argument has been discussed for mean curvature flows [36,10].Unfortunately, these results rely strongly on certain properties of the mean curvature flow (especially the comparison principle), which do not hold for fourth order equations like the Muskat problem.
There are also results in the literature [22,31] which eliminate the assumption for fourth order curvature driven flows (specifically the Stefan problem and the Mullins-Sekerka flow respectively).
However, the solutions constructed in [22] are discontinuous in time (the interface may experience sudden jumps in time) and the formulation in [31] only keeps track of the regular part of the interface.Let us also note that the Stefan problem and the Mullins-Sekerka flow are not similar to the Muskat problem.In particular, the jump condition for the pressure across the interface is different for the Muskat problem, which leads to qualitatively different behavior.

Paper summary
The rest of the paper is organized as follows.In Section 2, we derive the basic properties of the minimizing movements scheme (9) and construct our discrete-time quantities.We begin by showing that solutions to the minimization problem are characteristic functions at every time step of the discrete scheme.We then derive the existence of pressure as a Lagrange multiplier for the incompressibility constraint and obtain the Euler-Lagrange equation for the minimization problem.Our derivation of these equations were inspired by previous results from [11,25,19,18].In particular, the definition of the discrete in time pressure variable is very much inspired from [25].
In Section 3, we take τ and ε to zero together to obtain weak solutions to the Muskat problem, under the assumption that the internal energy of the discrete solutions converges to the internal energy of the limiting solutions.The main task in this section amounts to showing that that one can pass to the limit in the Euler-Lagrange equation obtained in Section 2. This can be done using the standard theory for Wasserstein gradient flows if ε is held fixed.However, the joint limit, τ, ε τ → 0, requires an adaptation of the arguments of [20] to the case of Wasserstein gradient flows.In particular, we develop a careful Aubin-Lions type argument, by applying a version of flow-exchange technique introduced in [24].Thus, while we largely follow ideas from [20] in our analysis, in some cases we give simplified proofs of the results we need, and in others we develop completely new arguments to conclude our results.
Finally, in Section 4, we conclude the paper with a discussion of the global minimizers of the approximated internal energy associated to the Muskat problem.While this discussion remains at the heuristic level, we conjecture the validity of the geometric description of the equilibrium configurations.
2. The Wasserstein minimizing movements scheme for the heat content 2.1.Some preliminary results.Recall that the setting for our problem is a smooth convex domain Ω ⊂ R d with L d (Ω) = 2.By P(Ω) we denote the space of Borel probability measures on R d supported on Ω. P ac (Ω) stands for the elements of P(Ω) that are absolutely continuous with respect to L d Ω.
Let Φ 1 , Φ 2 : Ω → R be given Lipschitz potentials and let us recall the definition of the potential energy Φ : P(Ω) × P(Ω) → R, given as Let ε > 0. We consider the heat content HC ε : P(Ω) × P(Ω) → R defied in (8), using the standard heat kernel We also use the notations We have the following preliminary results.Lemma 2.1.Let HC ε be defined as in (8).We have the following properties.
To show (2), it is enough to show that the function z → K ε (z) is λ-convex in the classical sense for some λ ∈ R. We have Since the matrix x ⊗ x is nonnegative definite for any We conclude similarly as in [11,Lemma 2.1] the displacement 2λ-convexity of HC ε on P(Ω) × P(Ω).
Proof.For any pair (ρ 1 , ρ 2 ) ∈ P we may write ρ 2 (x) = 1 − ρ 1 (x).Thus, extending the densities by 0 outside of Ω, we have Ignoring the constant multiples, both terms can be expressed conveniently in the Fourier domain: Now the strict concavity follows immediately as Ĝ(ξ 2.2.The minimizing movements scheme.Now we are ready to discuss the minimizing movements scheme.Our first result confirms the existence of minimizers, and shows that any minimizing configuration ρ = (ρ 1 , ρ 2 ) is a completely unmixed partition of the domain.As we will see, the phases stay unmixed thanks to the concavity of the heat content.
Proposition 2.3.Suppose that ρ n 1 , ρ n 2 ∈ P(Ω) and let τ > 0. Then the set of minimizers of the problem is non-empty and any solution (ρ * 1 , ρ * 2 ) ∈ P(Ω) × P(Ω) is the characteristic function of a partition of Ω. Proof.The existence of a solution of the optimization problem is an easy consequence of the weak lower semicontinuity of the objective functional and the (weak-⋆) compactness of P(Ω) × P(Ω).Let us remark that the constraint ρ 1 + ρ 2 = 1 a.e. is closed under weak convergence, since ´Ω ρ 1 dx + ´Ω ρ 2 dx = L d (Ω).
To show that an arbitrary solution (ρ * 1 , ρ * 2 ) is the characteristic functions of a partition of Ω, let us rewrite equivalently the minimization problem in terms of transport plans π i ∈ P(Ω × Ω).Recall that π i is a plan between ρ n i and ρ i , whenever ˆΩ×Ω ϕ(x) dπ i (x, y) = ˆΩ ϕ(x) dρ n i (x) and for any ϕ, ψ ∈ C(Ω).Moreover, since we are always working with measures ρ n i , ρ i that are absolutely continuous w.r.t.L d Ω, we can restrict our search to plans that have absolutely continuous marginals w.r.t.L d Ω.For such measures we will have for any ϕ, ψ ∈ C(Ω).For a measure θ ∈ P(Ω × Ω), we use the notation θ 1 := (P x ) # θ and θ 2 := (P y ) # θ to denote its marginals (here P x , P y : Ω × Ω → Ω stand for the canonical projections from Ω × Ω onto Ω).Thus, we aim to solve where we denote and we define and where we have extended the second marginals of π i by 0 outside of Ω.
The minimization is carried out over a weakly compact set and S is weakly lower semicontinuous and bounded below, thus minimizers exist.If π * = (π 1 , π 2 ) is a minimizer of ( 13) then we can construct a minimizer ρ * = (ρ * 1 , ρ * 2 ) of the original problem by taking ρ * i = (P y ) # π * i .Now we consider the properties of minimizers.Clearly, F is Gâteaux differentiable at π * in the sense that there exists δF where π + tθ is any admissible perturbation of π.Similarly, as the other terms in the definition of S are linear in π, these are in the same way differentiable, therefore S is Gâteaux differentiable in this sense.From Lemma 2.2 it follows that S is concave along line segements π + tθ (t ∈ (−1, 1)), where π is a feasible point and θ = (θ 1 , θ 2 ) is a feasible direction at π i.e. ( 15) and for some δ > 0 (16) in the sense of signed measures, for all t ∈ [0, δ).Furthermore, it follows from Lemma 2.2 that S is strictly concave on line segments π + tθ, if for some i the marginal θ 2 i (y) is not 0 for almost every y.Therefore, if π * = (π * 1 , π * 2 ) is a minimizer and θ is a feasible direction at π * with at least one non-trivial marginal then δS(π * ), θ > 0, where δS(π * ) stands for the first variation of S at π * defined as is (14).
Let π be a feasible solution and let ρ i (y) = π 2 i (y).Now suppose that there exists 0 < α < 1 such that the set Ω α = {y ∈ Ω : ρ 1 (y), ρ 2 (y) ∈ (α, 1 − α)} has positive measure.Partition Ω α into two sets E 1 , E 2 of equal measure.Then there exist measure preserving maps T : E 1 → E 2 and S : E 2 → E 1 such that S • T and T • S are the identity almost everywhere on their respective domains (for example one may choose T = ∇ψ to be the optimal transport map between the densities 1 E1 and 1 E2 and S = ∇ψ * ).Now we construct a feasible direction θ at π as follows.Let r(y) = ρ1(T (y)) ρ2(y) for y ∈ E 1 and we define the signed measures θ 1 , θ 2 ∈ M (Ω × Ω) as and for any ϕ ∈ C(Ω × Ω).Since θ 2 1 (y) = ρ 1 (T (y)) > α a.e. on E 1 we see that θ has a nontrivial marginal.Let us now check that θ is feasible, i.e. it satisfies (15) and ( 16).If β < min{1 − α, α} then π ± βθ defines a non-negative measure.Next, we check that θ satisfies θ 1 i (x) = 0 for a.e.x ∈ Ω and i = 1, 2 and where we have used that both T and S are measure preserving between E 1 and E 2 and vice-versa, respectively.Let us notice that these arguments also show (by taking ϕ ≡ 1) that Now, for ϕ ∈ C(Ω), we have Note that our arguments in fact show that −θ is also a feasible direction at π. ±θ have nontrivial marginals, and it is not possible to have both δS(π * ), θ > 0 and − δS(π * ), θ > 0. Therefore, π cannot be a minimizer.This allows us to conclude that for any minimizer ρ * of the original problem each density ρ * i takes values {0, 1} almost everywhere.

2.3.
Optimality conditions and construction of the pressure variables.In the next Lemma, we give a more complete characterization of the minimizers in terms of certain necessary inequalities.In particular, this is the first place where we see the appearance of the pressure variable, which plays an essential role in all of the subsequent analysis.Note that for convenience we express this result using the 2 ) be an optimizer in (11) and let (ρ 1 , ρ 2 ) a pair of probability measures such that ρ 1 + ρ 2 = 1 a.e.Then, (i) we have the following optimality condition for a suitable pair of Kantorovich potentials (ϕ 1 , ϕ 2 ) in the optimal transport of ρ * 1 onto ρ n 1 and ρ * 2 onto ρ n 2 , respectively.
(ii) there exists a Lipschitz continuous function p : Ω → R such that for any probability densities ρ 1 , ρ 2 with ρ 1 + ρ 2 = 1 a.e.we have Proof.Let (ρ 1 , ρ 2 ) be a pair of probability measures such that ), which by construction satisfy the constraint.
By optimality, we have Using the exact same argument as in [25,Lemma 3.1] to develop the Wasserstein part on the one hand and the first variations of HC ε and Φ on the other hand, we find (17).
For (ii) (similarly as in [19,Proposition 4.7]), let us notice first that ( 17) can be written in the form We know from Proposition 2.3 that ρ * 1 , ρ * 2 forms a partition of Ω.Therefore, we must take h 1 ≤ 0 on spt(ρ * 1 ), and h 1 ≥ 0 on spt(ρ * 2 ) (and vice-versa for h 2 ).To preserve the constraint ρ * 1 + δh 1 + ρ * 2 + δh 2 = 1 and the mass of each density, we must also take h 1 + h 2 = 0 a.e. and ´Ω h 1 dx = ´Ω h 2 dx = 0. Now if we set h 2 = −h 1 , we find that for any h 1 ∈ L ∞ (Ω) with 0 mean such that h 1 ≤ 0 a.e. on spt(ρ 1 ).This implies that there exist constants From here, we can define the pressure variable as (20) p e. on spt(ρ * 2 ).Since the Kantorovich potentials are Lipschitz continuous and the other terms are smooth, we find that p is Lipschitz continuous (note this regularity may degenerate as τ ↓ 0).By construction, p clearly satisfies the inequality in (18), and in particular the value of the l.h.s. is equal to zero.Since the functions under consideration are all Lipschitz continuous, we obtain (19).
2.4.Continuous in time solutions for ε > 0 fixed.Now we are ready to begin constructing time interpolations from the discrete scheme.In this section we restrict ourselves to the case where ε is held fixed.In this special case, we can use standard arguments from the theory of Wasserstein gradient flows to obtain continuous in time equations in the limit τ ↓ 0. When ε is held fixed, we have strong compactness on the pressure term.However, similarly to the models in [18,19], we will lack strong compactness on the density variables, which would mean that in the limit when τ ↓ 0, only a weaker version of the system will be available (see (26)).Let us also note that later on in Section 3, we will need some of the following pressure estimates when we take ε to zero along with τ .
Let T > 0 be a given time horizon and N ∈ N and τ > 0 such that N τ = T. Let ε > 0 be fixed.By now, it is standard how to construct weak solutions to PDEs that have gradient flow structure, using minimizing movement schemes as ( 21) Let ρ τ i : [0, T ] → P(Ω) be defined as ( 22) We define the corresponding velocities, pressures and momentum variables as ( 23) where ϕ n+1 i and p n+1 i are defined in Lemma 2.4.Following the same steps as in [18,25], the analysis boils down to obtain a sufficient amount of uniform (in τ ) estimates and compactness for the previously obtained functions, then pass to the limit τ ↓ 0. Also, as additional tools we construct the corresponding continuous in time (geodesic) interpolations between the densities and the corresponding velocities and momentum variables, (ρ τ i , ṽτ i , Ẽτ i ).We refer to [18, Section 3.1] (see also [25,33]) for the precise construction.In particular, as a consequence of this construction, we have Based on the same techniques as in [25,18], it is easy to obtain the following estimates.
and up to passing to subsequences they have the same distributional limits as τ ↓ 0.
(vi) Proof.Let us notice that the proofs of point (i), (ii) and (iii) follow the exact same lines and the proofs of [25, Lemma 3.3] and [18, Lemma 3.3], so we omit them.From Proposition 2.3 we know that that ρ n i 's are characteristic functions, therefore the proofs of (iv) and (v) follow the exact same lines as the proof of [20,Lemma 2.4.].
Let us give the details on (vi).From the identities (19) from Lemma 2.4, we have that there exists C > 0 (independent of τ > 0, which might increase from one inequality to the next) such that where we have used the uniform bounds on (ρ τ i ) τ >0 from the previous points and Taking the inner product of both sides in (19) with ξ(t, •) and integrating on Ω w.r.t.ρ n+1 i (we drop the dependence on t in the notation of ξ), we obtain ( 24) and interchanging the roles of ρ n+1 First, we have where in the second and third equalities we have used the change of variables x → x − √ εz and z → −z, respectively and the fundamental theorem of calculus in the last equality.Now, since for some universal constants α, β > 0, by the previous chain of equalities, we obtain where, in the last inequality we used the monotonicity of HC ε along the sequence (ρ n 1 , ρ n 2 ) n .
Furthermore, we have that Using the fact that we have piecewise constant interpolations, integrating the last equality in time on [0, T ], we get Now, adding together ( 24) and ( 25), then integrating in time on [0, T ] and using (ii), we obtain where the constant C > 0 is independent of τ > 0 and ε > 0. The thesis of follows.
Now we can use the above pressure estimates to derive a system of continuous in time equations in the limit τ → 0 when ε is held fixed.
Since p τ − ffl Ω p τ (•, x) dx τ >0 is uniformly bounded in L 2 ([0, T ]; H 1 (Ω)), there exists p ∈ L 2 ([0, T ]; H 1 (Ω)) such that up to passing to a subsequence, p τ − ffl Ω p τ (•, x) dx ⇀ p, weakly in L 2 ([0, T ]×Ω) and ∇p τ ⇀ ∇p, We only need to identify the limits of the momentum variables (E τ i ) τ >0 .From the weak convergence of the density variables, we have that up to passing to a subsequence, Last, by the previous arguments, we have Therefore, the thesis of the theorem follows.
It remains open whether in the previous theorem we have strong convergence ρ τ i → ρ i in L 2 ([0, T ] × Ω) as τ ↓ 0, and in particular whether one can claim that It is also open whether the continuum in time densities in Theorem 2.1 are characteristic functions.Indeed, since we can only guarantee that the discrete densities converge weakly to the continuum densities, the characteristic function property may be lost in the limit.We point out that due to the energy bounds, the densities are "almost" characteristic functions, i.e. we have Lemma 2.6.For (ρ τ 1 , ρ τ 2 ) given as above with ε > 0 fixed, for any α ∈ (0, 1) we have Proof.Let us set i = 1, the other case will be parallel.We have By Jensen's inequality, the above is smaller than We then estimate . Applying these inequalities, the result follows.

Muskat flow with surface tension
In this section, we complete the proof of Theorem 1.1 and show that when ε = ε τ goes to zero along with τ , the time interpolated minimizing movements scheme constructed in ( 22)-( 23) converges to a weak formulation of the Muskat Problem with surface tension under the energy convergence assumption (EC).This amounts to showing that each quantity in the system of Euler-Lagrange equations given in (19) converges (in an appropriate sense) to the correct limiting object.In particular, we will need strong L 1 convergence of the density functions in [0, T ] × Ω.To obtain the necessary compactness for strong L 1 convergence, we develop an adaptation of an Aubin-Lions type lemma in Proposition 3.2.We then conclude our result by verifying the convergence of the Euler-Lagrange equations to the weak formulation of the Muskat problem in a similar manner to the approach in [20].
Before we introduce the weak formulation of the Muskat Problem, let us recall the classical formulation of the problem.When the two phases are separated by a smooth interface Γ, the Muskat Problem is given by the continuity equation 0 along with the free boundary problem for the pressure where κ is the mean curvature of Γ, oriented to be positive when spt(ρ 2 ) is convex at the point.The weak formulation of the Muskat Problem with surface tension is provided in Definition 3.1 below.
We postpose the proof of the previous theorem to the end of this section.First, let us show a consistency result, i.e. that classical solutions of the Muskat problem satisfy the weak formulation ( 27) and ( 28).Lemma 3.1.If smooth solutions of (MP 1 )and (MP 2 ) exist with ρ i = χ Ai and with a C 2 hypersurface Γ = ∂A 1 = ∂A 2 , then ρ i satisfies ( 27) and (28) with the choice of (29) or using test functions Remark 3.2.Note that our notion of solution requires adding the surface measure σ 2 dµ to the classical pressure variable.This singular term appears from the minimizing movement scheme, where it ensures that a vacuum does not form at the interface.In general, we expect that the singular part in the weak pressure variable corresponds to the surface measure in (29).
Proof. ( 27) is a standard weak expression of (MP 1 ) with b i v i = −∇p i − ∇Φ i .Let us write again Γ = ∪ t>0 (Γ t × {t}).Then we have Here for the second equality we used integration by parts for the first integral, using the fact that and for the third equality we used the curvature jump condition at the interface.To conclude, let us recall that κ is oriented to be positive when spt(ρ 2 ) is convex at the point.With this orientation, observe that for any C 1 vector field ξ we have where ν = Dρ 1 /|Dρ 1 | is normal to Γ toward the support of ρ 1 .Hence we conclude that (28) holds.
Let us notice that in order to use the Aubin-Lions argument, we need to show a tightness condition (cf.[32, Definition 1.3, Remark 1.5]) for the time-dependent family (ρ τ i ) τ >0 .This is typically done by providing a compact set (of the space of probability measures), where 'most' of the sequences (ρ τ i (t)) τ >0 lie.Note that the estimate in Lemma 2.5(v) does not provide such a compact set.Indeed, for τ > 0, the densities (ρ τ i (t)) τ >0 are not actually BV in space.Inspired by arguments in [12], we will work first with an auxiliary sequence where we have performed a convolution with the heat kernel G ε only in the spacial variable.It worth noticing that this 'perturbation' of ρ τ i is reminiscent to the one obtained via the so-called flow interchange technique introduced in [24].We have the following properties for this new sequence.

Claim 1. (ρ τ
i ) τ >0 is uniformly (w.r.t.ε and τ ) in L 1 ([0, T ]; BV (Ω)).Proof of Claim 1.The uniform L ∞ bounds on (ρ τ i ) τ >0 are also preserved for (ρ τ i ) τ >0 .Now, as in the proof of [20,Lemma 2.4] we conclude that for a uniform constant C, which proves our claim.Claim 2. There exists a subsequence of (ρ τ i ) τ >0 (that we do not relabel) and For t ∈ (0, T ) fixed, let us notice that ρ τ i,τ actually can be seen as the evolution of ρ τ i,t via the heat flow for time ε > 0. It is well-known that W 2 is contractive along the heat flow, so we have By construction, F is convex, l.s.c. in L 1 (Ω) and it sublevel sets are compact in L 1 (Ω), therefore it defines a normal coercive integrand.
From (30) and from the definition of F , we have And similarly from Lemma 2.5(i) and from (31), we have lim therefore the assumptions of [32,Theorem 2] are fulfilled and one can conclude that there exists ρ i ∈ L 1 ([0, T ] × Ω) and a subsequence of (ρ τ i ) τ >0 (that we do not relabel) which is converging in measure, and in particular pointwise a.e. to ρ i as ε τ ↓ 0. The strong convergence in L 1 ([0, T ] × Ω) follows from Lebesgue's dominated convergence theorem, since (ρ τ i ) τ >0 is uniformly bounded.The claim follows.Claim 3.There exists a subsequence of the original sequence (ρ τ i ) τ >0 which is converging to ρ i strongly in L 1 ([0, T ] × Ω) as ε τ ↓ 0.
Proof of Claim 3. Passing to the same subsequence in the original sequence (ρ τ i ) τ >0 (that we do not relabel) as in the last step of the proof of Claim 2, we have , for a constant C > 0 (independent of τ and ε) and the claim follows by taking ε τ ↓ 0. In the last inequality we have used where we relied on the fact that since ρ τ 1 and ρ τ 2 take values in {0, 1}, we have To conclude, let us notice that since (ρ τ i ) τ >0 converges uniformly w.r.t.W 2 to ρ i , ρ i and ρ i have to coincide.Also, for this last subsequence, we can pass to the limit as ε τ ↓ 0 in the estimation of Lemma 2.5(v) to obtain that actually ρ i ∈ L 1 ([0, T ]; BV (Ω; {0, 1})).
To overcome this issue, let us argue using the following claim.This is a consequence of classical result in the theory of elliptic equations and Schauder estimates (see for instance [ where in the first inequality we have used the last uniform estimate from the proof of Lemma 2.5(vii).This implies that (p τ (t, •)) τ is uniformly bounded in (C 0,α (Ω)) * .
To conclude the thesis of the lemma, we observe that (by possibly subtracting the mean) (p τ ) τ is also converging weakly-⋆ to some p, therefore we will have that ζ = ∇p in D ′ ((0, T ) × Ω).
To complete the proof of of Theorem 1.1, it remains to show that limit of equation (32) converges to the weak formulation of the Muskat problem.This result is proven in Proposition 3.3, which in turn uses the results of Lemmas 3.4 and 3.5.Note that the energy convergence assumption is finally used in Lemma 3.4.The assumption plays a crucial role in the following arguments as it allows us to convert a limit requiring uniform convergence to one which only requires pointwise convergence.In what follows, it will be useful for notational simplicity to introduce the following definition: Definition 3.2.For a smooth rapidly decaying kernel J : R d → R and a vector valued Radon measure ν ∈ M d (Ω) we define the Radon measure σ J (ν) ∈ M (Ω) as Proposition 3.3.Let (ρ τ i ) τ >0 be the piecewise constant interpolations constructed in ( 22)-( 23) and let ρ i ∈ L 1 ([0, T ]; BV (Ω; {0, 1})), i = 1, 2 be their strong L 1 limits.If the assumption (EC) is fulfilled, then up to passing to a subsequence that we do not relabel, we have Here, we denoted J(z) := σ √ 2π|z • e 1 | 2 G(z) and we used the decomposition ∇ξ sym is an appropriate basis of the the space of symmetric matrices.
as ε τ ↓ 0. This result is not straight forward, since both the functional and the densities depend on τ .This difficulty is the key reason that we need the energy convergence assumption (EC).Arguing exactly as in the proof of [20,Lemma 2.8], in order to show the convergence result (33), it is enough to show its time-independent version, i.e.
A computation similar to the one in the proof of Lemma 2.5(viii) reveals where, we were using a second order Taylor expansion (and the smoothness of ξ) in the last equality.Let us observe that the very last term in ( 34) is converging to 0 as ε τ ↓ 0 (due to the strong convergence ρ τ i → ρ i in L 1 (Ω) and the fact that ρ 1 ρ 2 ≡ 0 a.e., see Proposition 3.2).Therefore, it remains to prove as ε τ ↓ 0. We note that where e i ∈ R d is the i th standard basis vector.All of these basis matrices have the form n ⊗ n for some n ∈ S d−1 .Thus, we may write where {n k } is any indexing of the above matrices and ζ k (x)n k ⊗ n k = P k ∇ξ sym (x) where P k is the appropriate change of basis matrix.Therefore we have σ as ε τ → 0, and so, in this case we would have This would in fact complete the claim, as we can compute where we have used Tr(n k ⊗ n k ) = 1 and the fact that the antisymmetric parts of ∇ξ are annihilated by the above operations.
The conclusion in the previous proposition is made by the following two lemmas.where we define the Radon measure σ J (Dρ i ) ∈ M (Ω) as in Definition 3.2 and we have used the notation J ε (z) := J(z/ε).
The proof supplied below is different from the one by Laux-Otto in [20, Lemma 2.8, Lemma 3.6].Instead of disintegrating on R d and using one-dimensional arguments we obtain upper and lower bounds using mollifiers.
Proof of Theorem 3.1.This proof is a direct consequence of Theorem 3.2 and Proposition 3.3.
where Φ i (x) = c i x d , with 0 < c 1 < c 2 .The order of the constants denote that ρ 1 is the lighter fluid, where the vector −e d denotes the direction of gravity.The coordinate here is x = (x ′ , x d ) and for simplicity we consider a cylindrical domain, Here the convolution is taken with the extension of the density function as zero outside of the domain, with the density constraint ρ 1 + ρ 2 = 1, ´Rd dρ i = M i .
Let us mention that away from the global equilibrium, there are diverse possibilities of stationary states for Ẽε even with zero potentials.For instance any choice of characteristic functions ρ 1 and ρ 2 generating the interface as a disjoint union of spheres, {|x − a i | = r i }, is a stationary solution of the limit energy.
In the limit ε → 0, the Γ-convergence properties indicate that the global equilibrium of the ε-energy converges to the limiting density pair (ρ 1 , ρ 2 ) = (χ A , χ A c ), which is the global minimizer of the limit energy When ∂A is away from the lateral and bottom portion of the cylinder, this corresponds to the classical pendant liquid droplet problem where the minimizer is known to be rotationally symmetric and convex (see [15]).When the droplet boundary touches the cylinder boundary, various shapes of drops are possible (see Figure 1) and the complete description of possibly non-smooth global minimizers appear to be open.In general, when Σ has C 1,α boundary, it is shown in [37] that ∂A is C 1,α up to the boundary and meets ∂Σ orthogonally.For further discussion of available results we refer to [23].An ongoing work on numerical simulations for our flow suggest that the equilibrium states even for the ε-energy can be categorized as the ones on Figure 1.In dimension two, we could observe an additional equilibrium shape, as in Figure 2.

Figure 1 .
Figure 1.Equilibrium shapes in any dimensions, depending on the proportion of the lighter vs. the heavier fluid: orange stands for the lighter fluid, while gray stands for the heavier fluid

Figure 2 .
Figure 2.An additional equilibrium shape in dimension two are rotation invariant, when integrating on R d , it is enough to consider only n k = e 1 , where e 1 is the first element of the standard basis on R d .Now, defining J(z) := σ √ 2π|z • e 1 | 2 G(z), we deduce the claim in (iii) from Lemma 3.4 and Lemma 3.5.