Weak solutions to the Muskat problem with surface tension via optimal transport

TL;DR

This paper introduces a novel optimal transport-based framework to approximate weak solutions of the Muskat problem with surface tension, enabling convergence analysis and numerical simulations.

Contribution

It presents a new gradient flow approach in Wasserstein space for the Muskat problem, using heat content energy to relax surface tension effects and prove convergence.

Findings

Convergence of the scheme to weak solutions under energy assumptions

Numerical simulations demonstrating the scheme's effectiveness

Analysis of equilibrium configurations

Abstract

Inspired by recent works on the threshold dynamics scheme for multi-phase mean curvature flow (by Esedo\={g}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\={g}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.

Figures

Equations

$$v_i+b_i^{-1}\nabla {\delta }_{{\rho }_i}\mathcal{E}(\mathbf{\rho })=0,$$

$${\partial }_t{\rho }_i+\nabla \cdot ({\rho }_i{v}_i)=0,$$

$$\mathcal{E}(\mathbf{\rho })={\mathcal{E}}_p(\mathbf{\rho })+{\mathcal{E}}_s(\mathbf{\rho })+\Phi (\mathbf{\rho }).$$

$${\mathcal{E}}_p(\mathbf{\rho })=\{\begin{array}{l}0,\text{if}{\rho }_1(x)+{\rho }_2(x)=1\text{for a.e.}x\in \Omega \\ +\infty ,\text{otherwise}.\end{array}$$

$${\mathcal{E}}_s(\mathbf{\rho })=\{\begin{array}{l}\frac{\sigma }{2}|D{\rho }_1|(\Omega )+\frac{\sigma }{2}|D{\rho }_2|(\Omega ),\text{if}{\rho }_1,{\rho }_2\in BV(\Omega ;\{0,1\})\mathrm{and}{\rho }_1(x){\rho }_2(x)=0\text{for a.e.}x\in \Omega \\ +\infty ,\text{otherwise},\end{array}$$

$$\Phi (\mathbf{\rho })={\int }_{\Omega }{\Phi }_1\mathrm{d}{\rho }_1+{\int }_{\Omega }{\Phi }_2\mathrm{d}{\rho }_2,$$

$${\Phi }_i=g_{i}x\cdot e_d,g_i>0,i=1,2$$

$$\partial_{t}\rho_{i}-b^{-1}_{i}\nabla\cdot\big{(}(\nabla p_{i}+\nabla\Phi_{i})\rho_{i}\big{)}=0$$

$$\left\{\begin{array}[]{lll}-\Delta p_{i}=\Delta\Phi_{i}&\hbox{ in }&\operatorname{spt}(\rho_{i});\\ \\ \partial_{n}(p_{i}+\Phi_{i})=0&\hbox{ on }&\partial\Omega;\\ \\ V=b_{1}^{-1}\partial_{n}(p_{1}+\Phi_{1})=b_{2}^{-1}\partial_{n}(p_{2}+\Phi_{2})&\hbox{ on }&\Gamma;\\ \\ \left[p\right]:=(p_{1}-p_{2})=\frac{\sigma}{2}\kappa&\hbox{ on }&\Gamma;\\ \\ \tilde{n}=n&\hbox{on}&\partial\Gamma\cap\partial\Omega;\\ \end{array}\right.$$

$${\mathbf{\rho }}^{n+1}=\underset{\mathbf{\rho }}{\mathrm{argmin}}\left\{\mathcal{E}(\mathbf{\rho })+\sum _{i=1}^2\frac{b_i}{2\tau }{W}_2^2({\rho }_i,{\rho }_i^n)\right\}$$

$${\text{HC}}_{\epsilon }(\mathbf{\rho }):=\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{\Omega }(G_{\epsilon }\star {\rho }_1)(x)\mathrm{d}{\rho }_2(x)=\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{\Omega }{\int }_{{\mathbb{R}}^d}G(z)\mathrm{d}{\rho }_1(x+\sqrt{\epsilon }z)\mathrm{d}{\rho }_2(x).$$

$${\mathbf{\rho }}^{n+1}=\underset{\mathbf{\rho }}{\mathrm{argmin}}\left\{{\mathcal{E}}_{\epsilon }(\mathbf{\rho })+\sum _{i=1}^2\frac{b_i}{2\tau }{W}_2^2({\rho }_i,{\rho }_i^n)\right\},$$

$${\mathcal{E}}_{\epsilon }(\mathbf{\rho }):={\mathcal{E}}_p(\mathbf{\rho })+{\mathrm{HC}}_{\epsilon }(\mathbf{\rho })+\Phi (\mathbf{\rho }).$$

$$\frac{1}{N^2}\sum _{i\in P_1,j\in P_2}V(|x_i-x_j|)=\frac{1}{N^2}\sum _{i\in P_1,j\in P_2}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}V(|x-x^{\prime }|)\delta (x-x_i)\delta (x^{\prime }-x_j)$$

$${\rho }_1(t,\cdot )={\chi }_{A_t\cap \Omega }and{\rho }_2(t,\cdot )={\chi }_{\Omega \setminus A_t},$$

$${\int }_0^T{\mathrm{HC}}_{\epsilon }({\rho }_1^{\epsilon ,\tau },{\rho }_2^{\epsilon ,\tau })\mathrm{d}t\to {\int }_0^T\left(\sum _i\frac{\sigma }{2}{\int }_{\Omega }|D{\rho }_i|\right)\mathrm{d}x\mathrm{d}t,$$

$$\Phi (\mathbf{\rho }):={\int }_{\Omega }{\Phi }_1\mathrm{d}{\rho }_1+{\int }_{\Omega }{\Phi }_2\mathrm{d}{\rho }_2.$$

$$G_{\epsilon }(x)=\frac{1}{(4\pi \epsilon {)}^{d/2}}{e}^{-\frac{|x{|}^2}{4\epsilon }}.$$

$$D^2{K}_{\epsilon }(x)=\frac{\sigma \sqrt{2\pi }}{2(4\pi {)}^{d/2}}\frac{1}{{\epsilon }^{(d+3)/2}}{e}^{-\frac{|x{|}^2}{4\epsilon }}\left(\frac{1}{2\epsilon }x\otimes x-I_d\right).$$

$${\mathrm{HC}}_{\epsilon }({\rho }_1,{\rho }_2)=\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{{\mathbb{R}}^d}(G_{\epsilon }\star {\rho }_1)(x)(1-{\rho }_1(x))\mathrm{d}x=\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{{\mathbb{R}}^d}(G_{\epsilon }\star {\rho }_1)(x)\mathrm{d}x-\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{{\mathbb{R}}^d}(G_{\epsilon }\star {\rho }_1)(x){\rho }_1(x)\mathrm{d}x$$

$${\hat{\rho }}_1(0)-{\int }_{{\mathbb{R}}^d}\hat{G}(\xi \sqrt{\epsilon })|{\hat{\rho }}_1(\xi ){|}^2\mathrm{d}\xi .$$

$$\inf \left\{{\mathrm{HC}}_{\epsilon }(\mathbf{\rho })+\Phi (\mathbf{\rho })+\frac{b_1}{2\tau }{W}_2^2({\rho }_1,{\rho }_1^n)+\frac{b_2}{2\tau }{W}_2^2({\rho }_2,{\rho }_2^n):{\rho }_1,{\rho }_2\in {\mathcal{P}}^{\mathrm{ac}}(\Omega ),{\rho }_1+{\rho }_2=1a.e.\right\}$$

$${\int }_{\Omega \times \Omega }\phi (x)\mathrm{d}{\pi }_i(x,y)={\int }_{\Omega }\phi (x)\mathrm{d}{\rho }_i^n(x)\mathrm{and}{\int }_{\Omega \times \Omega }\psi (y)\mathrm{d}{\pi }_i(x,y)={\int }_{\Omega }\psi (y)\mathrm{d}{\rho }_i(y),$$

$$\begin{array}[]{l}\displaystyle\inf\left\{\mathcal{F}(\pi_{1},\pi_{2})+\mathcal{G}(\pi_{1},\pi_{2})+\sum_{i=1}^{2}\int_{\Omega\times\Omega}\frac{b_{i}}{2\tau}|x-y|^{2}\,{\rm d}\pi_{i}(x,y)\right\}=:\inf\mathcal{S}(\pi_{1},\pi_{2})\\[5.0pt] \displaystyle{\rm{subject\ to\ }}\pi_{i}^{1}=\rho^{n}_{i}\ {\rm{and}}\ \sum_{i=1}^{2}\pi_{i}^{2}=1\ {\rm a.e.}\end{array}$$

$$\mathcal{G}({\pi }_1,{\pi }_2):={\int }_{\Omega \times \Omega }\sum _{i=1}^2{\Phi }_i(y)\mathrm{d}{\pi }_i(x,y).$$

$$\mathcal{F}({\pi }_1,{\pi }_2)$$

$$-\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{\Omega \times \Omega }{\int }_{\Omega \times {\mathbb{R}}^d}G(y_1-\sqrt{\epsilon }{y}_2)\mathrm{d}{\pi }_1(x_2,y_2)\mathrm{d}{\pi }_1(x_1,y_1)$$

$$\mathcal{S}({\pi }_1,{\pi }_2):=\mathcal{F}({\pi }_1,{\pi }_2)+\mathcal{G}({\pi }_1,{\pi }_2)+\sum _{i=1}^2{\int }_{\Omega \times \Omega }\frac{b_i}{2\tau }|x-y{|}^2\mathrm{d}{\pi }_i(x,y),$$

$$⟨\delta \mathcal{F}({\mathbf{\pi }}^{\ast }),\mathbf{\theta }⟩$$

$$-\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{\Omega \times \Omega }{\int }_{\Omega \times {\mathbb{R}}^d}G(y_1-\sqrt{\epsilon }{y}_2)\mathrm{d}{\pi }_1(x_2,y_2)\mathrm{d}{\theta }_1(x_1,y_1)$$

$$-\sigma \sqrt{\frac{2\pi }{\epsilon }}{\int }_{\Omega \times \Omega }{\int }_{\Omega \times {\mathbb{R}}^d}G(y_1-\sqrt{\epsilon }{y}_2)\mathrm{d}{\theta }_1(x_2,y_2)\mathrm{d}{\pi }_1(x_1,y_1),$$

$${\theta }_i^1(x)=0a.e.x\in \Omega ,{\int }_{\Omega \times \Omega }\mathrm{d}{\theta }_i(x,y)=0\mathrm{and}\sum _{i=1}^2{\theta }_i^2(y)=0a.e.y\in \Omega ,$$

Full text

1. Introduction

The Muskat problem was first introduced by Morris Muskat [30] as a model for the flow of two immiscible fluids through a porous medium. Since its introduction, this problem has received sustained attention in a variety of fields. It 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

[TABLE]

coupled with the continuity equation

[TABLE]

where $v_{i}$ is the velocity of phase $i$, $\bm{\rho}=(\rho_{1},\rho_{2})$ is the collection of relative concentrations for each phase, $\nabla(\delta_{\rho_{i}}\mathcal{E})$ denotes the spacial gradient of the classical first variation of the free energy with respect to $\rho_{i}$, and $b_{i}>0$ ($i=1,2$) denotes constant mobilities. For convenience, throughout the rest of the paper, we will refer to $\bm{\rho}$ as a collection of density functions; however, one should note that $\bm{\rho}$ 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 $\Omega\subset\mathbb{R}^{d}$ with smooth boundary. We shall suppose that the two fluids fill the entire domain, and that they are confined to $\Omega$ for all time. We then take the internal energy to be a sum of three distinct terms:

[TABLE]

The first term in the energy, describing incompressibility and containment of the fluids, is given by

[TABLE]

The immiscibility of the fluids and the surface tension force arise from the highly non-convex interaction energy

[TABLE]

where $|D\rho_{i}|(\Omega)$ denotes the total variation of $\rho_{i}$ in $\Omega$ and $\sigma>0$ is a surface tension constant. Finally, $\Phi(\bm{\rho})$ denotes the potential energy of the fluid configuration, i.e.

[TABLE]

where $\Phi_{1},\Phi_{2}:\Omega\to\mathbb{R}$ are given Lipschitz continuous potentials. A typical example is when one assumes these to be gravitational potentials, i.e.

[TABLE]

where $g_{i}$’s are proportional to the specific gravity of each fluid.

Although the internal energy is singular, when $\rho_{1}$ and $\rho_{2}$ are separated by a smooth interface $\Gamma:=\partial\{\rho_{1}>0\}\cap\partial\{\rho_{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 $\mathcal{E}_{p}$ above. The continuity equation becomes

[TABLE]

and the pressure is determined by solving the free boundary problem

[TABLE]

where $\kappa$ denotes the mean curvature of $\Gamma$, oriented to be positive when $\{\rho_{2}>0\}$ is convex at the point, $n$ denotes the outer normal along $\partial\Omega$ and along $\Gamma$, and $\tilde{n}$ denotes the co-normal vector orthogonal to $\partial\Gamma$ and tangential to $\Gamma$. Note that the final condition relating the co-normal vector at $\partial\Gamma\cap\partial\Omega$ to the normal vector of $\partial\Omega$ implies that $\Gamma$ must meet $\partial\Omega$ orthogonally (see Lemma 3.1 and Remark 3.2 for the weak formulation of this condition). To summarize the ideas in the formal derivation of (MP1)-(MP2) from (1)-(2) using the definition of $\mathcal{E}$, we heuristically have $\nabla\delta_{\rho_{i}}\mathcal{E}_{p}=\nabla p_{i}$, $\nabla\delta_{\rho_{i}}\Phi=\nabla\Phi_{i}$ while the contribution of $\nabla\delta_{\rho_{i}}\mathcal{E}_{s}$ will act only on $\Gamma$ in the form of the curvature $\kappa$.

Problems like (MP2) received a lot of attention in the past decades. Most of the works focus on the zero surface tension model ($\sigma=0$) and well-posedness of regular solutions with graph property [2, 6, 8, 9, 10]. In the presence of surface tension, the problem has stronger regularity properties in stable settings [34], but still, topological singularities can occur in finite time, for instance when heavier fluid is placed on top of the lighter one [14]. Thus, our aim is to construct global-in-time weak solutions to the Muskat problem (MP2), 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 [31, 32], Darcy’s law can be approximated by the Euler-Lagrange equation for the minimizing movements scheme (or JKO [19] scheme) with time step size $\tau>0$,

[TABLE]

where $W_{2}(\rho_{i},\rho_{i}^{n})$ 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. [26, 5]), 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 $\tau>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 non-convex constraint set $\{\rho_{1},\rho_{2}\in BV(\Omega;\{0,1\}):\rho_{1}(x)\rho_{2}(x)=0\;\textrm{a.e.}\}$ is extremely difficult. For these reasons, we instead consider a relaxed version of minimizing movements scheme inspired by [13] and [22].

Approximation of the perimeter by the Heat Content

In our analysis, we replace the interaction energy $\mathcal{E}_{s}$ in (7) by the heat content energy

[TABLE]

Here $G_{\varepsilon}:\mathbb{R}^{d}\to\mathbb{R}$ stands for the standard heat kernel (with mean 0 and variance $\varepsilon>0$) and the densities $\rho_{i}$ are assumed to be defined on all of $\mathbb{R}^{d}$ by extending them to zero off of $\Omega$. Let us notice that in [13] and [22], for similar purposes the authors use periodic extensions. The analysis in both cases and the validity of results using both kinds of extensions is essentially the same.

Dating back to the work of De Giorgi, approximate perimeter energies have been used in the literature to study geometric variational problems (see for instance [29] and [1]). The use of the heat content energy to study the multi-phase mean curvature flow was first introduced by Esedoḡlu and Otto in [13]. It was observed in [13] that the threshold dynamics, a well known numerical scheme for mean curvature motion introduced by Merriman, Bence, and Osher [28], is precisely a minimizing movements scheme for the heat content energy. Esedoḡlu and Otto also showed that $\textrm{HC}_{\varepsilon}$ $\Gamma$-converges (with respect to the $L^{1}$ topology) to $\mathcal{E}_{s}$ as $\varepsilon\to 0$. Building off of these results, Laux and Otto showed in [22] that under an energy convergence assumption, the threshold dynamics scheme produces weak solutions to the multi-phase motion by mean curvature in the limit $\varepsilon\to 0$.

Our goal is to consider such a framework in the context of the Muskat problem by studying the minimizing movements scheme

[TABLE]

where we used the notation

[TABLE]

As we alluded above, the scheme (9) has a number of numerical advantages over (7). Unlike $\mathcal{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 Subsection 5.1 of [13]. After applying this trick, the resulting variational problem becomes convex, and thus, can be efficiently solved using the recently introduced back-and-forth method [18]. See Figures 1-3, for a demonstration of the numerical performance of the scheme.

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 ${\rm{HC}}_{\varepsilon}$ 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. [17]). Typically, the discrete surface tension takes the form

[TABLE]

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$ ([17]). By taking the limit $N\to\infty$ in above formula, one obtains an analogue of the heat content energy where the kernel $G_{\varepsilon}$ 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 (just as in the case of [13] and [22]). This includes scenarios where the surface tension force depends on the phases that are interacting [13]. To present our ideas in the simplest possible way, we do not pursue the multiphase case in this paper.

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 $\tau$ and the heat content approximation parameter $\varepsilon$ 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. In addition, we show that our weak formulation (see Definition 3.1) encodes all of the conditions in (MP1)-(MP2) (see Lemma 3.1 and Remark 3.2).

Theorem 1.1**.**

Let $T>0$ be a fixed time horizon, let $\varepsilon,\tau>0$ be fixed and let $(\rho_{1}^{\varepsilon,\tau},\rho_{2}^{\varepsilon,\tau})$ be the discrete in time interpolations of the densities obtained from the minimizing movements scheme (9). Let moreover $p^{\varepsilon,\tau}$ stand for the discrete in time interpolations between the scalar pressure fields, obtained as Lagrange multipliers associated to the incompressibility constraint in (9). Then

• •

There exists a family $(A^{\varepsilon,\tau}_{t})_{t\in[0,T]}\subseteq\Omega$ of measurable sets such that $\rho_{1}^{\varepsilon,\tau}(t,\cdot)=\chi_{A^{\varepsilon,\tau}_{t}}$ and $\rho_{2}^{\varepsilon,\tau}(t,\cdot)=\chi_{\Omega\setminus A^{\varepsilon,\tau}_{t}}$.

• •

There exists $\rho_{i}\in L^{1}([0,T];BV(\Omega;\{0,1\}))\cap{\rm{AC}}^{2}([0,T];{\mathscr{P}}(\Omega))$ ($i=1,2$) such that as $\max\{\varepsilon,\tau\}\downarrow 0$ and along a subsequence $(\rho_{1}^{\varepsilon,\tau},\rho_{2}^{\varepsilon,\tau})\to(\rho_{1},\rho_{2})$ strongly in $L^{1}([0,T]\times\Omega)\times L^{1}([0,T]\times\Omega)$. Moreover, $\rho_{1},\rho_{2}\in L^{1}([0,T];BV(\Omega))$ and they are also characteristic functions that sum up to one, i.e.

[TABLE]

for a measurable family of sets $(A_{t})_{t\in[0,T]}$ which are of finite perimeter.

• •

There exists a scalar pressure field $p\in L^{2}([0,T];(C^{0,\alpha}(\Omega))^{*})$ such that along a subsequence $\nabla p^{\varepsilon,\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}\nabla p$ weakly-$\star$* in $L^{2}([0,T];(C^{1}(\Omega))^{*})$ as $\max\{\varepsilon,\tau\}\downarrow 0$.*

• •

Finally, under the assumption of energy convergence

[TABLE]

$(\rho_{i},v_{i},p)$* with $i=1,2$ solves, in the weak sense (see Definition 3.1), the problem (MP1)-(MP2). Here formally $v_{i}$ corresponds to $-b_{i}^{-1}(\nabla p+\nabla\Phi_{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. This is verified with several numerical simulations in Section 4.1.

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 $\Omega$. 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 $(\rho_{1}^{\varepsilon,\tau},\rho_{2}^{\varepsilon,\tau})$ in the case when $\varepsilon>0$ is fixed and $\tau\downarrow 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\{\tau,\varepsilon\}\downarrow 0$ we recover total variation estimates on the densities. This compactness is obtained via a standard Aubin-Lions type argument.

• (3)

The energy convergence assumption (EC) is rather natural and the same as the ones given in [22] and [23]. This assumption ensures that there is no sudden loss of boundary between phases in the limit $\varepsilon\to 0$. If there is a loss of interface then one cannot obtain weak solutions. Indeed, if two components of the support of $\rho_{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 $\Omega$, creating an inconsistency between the discrete and limiting evolutions.

Replacing the assumption with a direct argument has been discussed for mean curvature flows [40, 11]. 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 [24, 35] which eliminate the assumption for third order curvature driven flows (specifically the Stefan problem and the Mullins-Sekerka flow respectively). However, the solutions constructed in [24] are discontinuous in time (the interface may experience sudden jumps in time) and the formulation in [35] 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 [12, 27, 21, 20]. In particular, the definition of the discrete in time pressure variable is very much inspired by [27].

In Section 3, we take $\tau$ and $\varepsilon$ 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 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 $\varepsilon$ is held fixed. However, the joint limit, $\tau,\varepsilon_{\tau}\to 0$, requires an adaptation of the arguments of [22] to the case of Wasserstein gradient flows. Let us underline that one can rely entirely on the results of [22] to pass to the limit the weak curvature equation. An adaptation of the Aubin-Lions type argument from [22] can be done to get compactness of the density terms. The only difference here is that we are using $W_{2}$ as metric while in [22] a different metric is used, but this does not impose crucial difficulties. An interesting link to the flow-exchange technique introduced in [26], which is a typical tool for $W_{2}$ gradient flows, is also pointed out. Last, we developed the necessary estimates and compactness results on the pressure terms. These are new and clearly were not present in the setting of mean curvature flows. The compactness that we get on the pressure is in the sense of distributions.

Finally, in Section 4, we conclude the paper with a demonstration of the numerical method on several examples and a discussion on the global minimizers of the approximated internal energy associated to the Muskat problem. While this discussion remains at the heuristic level, our conjectures are supported by the equilibrium states attained in our numerical experiments (c.f. Figures 1-3). We end the paper with an appendix section, where we recall the results from [22] that are used when passing to the limit the weak curvature equation.

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 $\Omega\subset\mathbb{R}^{d}$ and without loss of generality by scaling we assume that $\mathscr{L}^{d}(\Omega)=2$. By ${\mathscr{P}}(\Omega)$ we denote the space of Borel probability measures on $\mathbb{R}^{d}$ supported on $\overline{\Omega}$. ${\mathscr{P}}^{{\rm ac}}(\Omega)$ stands for the elements of ${\mathscr{P}}(\Omega)$ that are absolutely continuous with respect to $\mathscr{L}^{d}\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}\Omega$.

Let $\Phi_{1},\Phi_{2}:\Omega\to\mathbb{R}$ be given Lipschitz potentials and let us recall the definition of the potential energy $\Phi:{\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega)\to\mathbb{R}$, given as

[TABLE]

Let $\varepsilon>0$. We consider the heat content ${\rm{HC}}_{\varepsilon}:{\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega)\to\mathbb{R}$ defied in (8), using the standard heat kernel $G_{\varepsilon}:\mathbb{R}^{d}\to\mathbb{R}$, i.e.

[TABLE]

We also use the notations $K_{\varepsilon},G:\mathbb{R}^{d}\to\mathbb{R}$ to denote $K_{\varepsilon}(x)=\sigma\sqrt{\frac{2\pi}{\varepsilon}}G_{\varepsilon}(x)$ and $G(x)=G_{1}(x).$ We have the following preliminary results.

Lemma 2.1**.**

Let ${\rm{HC}}_{\varepsilon}$ be defined as in (8). We have the following properties.

• (1)

${\rm{HC}}_{\varepsilon}$* is bounded from below and continuous w.r.t. the weak-$\star$ convergence on ${\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega).$*

• (2)

${\rm{HC}}_{\varepsilon}$* displacement $\lambda$-convex on ${\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega)$, with $\lambda=-\frac{\sigma\sqrt{2\pi}}{(4\pi)^{d/2}}\frac{1}{\varepsilon^{(d+3)/2}}$.*

Proof.

(1) Is immediate by the definition of ${\rm{HC}}_{\varepsilon}$.

To show (2), it is enough to show that the function $z\mapsto K_{\varepsilon}(z)$ is $\lambda$-convex in the classical sense for some $\lambda\in\mathbb{R}$. We have

[TABLE]

Since the matrix $x\otimes x$ is positive semidefinite for any $x\in\mathbb{R}^{d}$, setting $\lambda=-\frac{1}{2(4\pi)^{d/2}}\frac{\sigma\sqrt{2\pi}}{\varepsilon^{(d+3)/2}}$, we have that the matrix $D^{2}K_{\varepsilon}(x)-\lambda I_{d}$ is positive semidefinite for any $x\in\mathbb{R}^{d},$ which implies in particular that $K_{\varepsilon}$ is $\lambda$-convex.

We conclude similarly as in [12, Lemma 2.1] the displacement $2\lambda$-convexity of ${\rm{HC}}_{\varepsilon}$ on ${\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega).$ ∎

Lemma 2.2**.**

If $P\subset{\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega)$ denotes the set of pairs $(\rho_{1},\rho_{2})$ such that $\rho_{1}+\rho_{2}=1$ a.e. in $\Omega$, then the heat content is strictly concave along line segments in $P$.

Proof.

For any pair $(\rho_{1},\rho_{2})\in P$ we may write $\rho_{2}(x)=1-\rho_{1}(x)$. Thus, extending the densities by [math] outside of $\Omega$, we have

[TABLE]

Ignoring the constant multiples, both terms can be expressed conveniently in the Fourier domain:

[TABLE]

Now the strict concavity follows immediately as $\hat{G}(\xi\sqrt{\varepsilon})>0$ for all $\xi\in\mathbb{R}^{d}$. ∎

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 $\bm{\rho}=(\rho_{1},\rho_{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 $\rho_{1}^{n},\rho_{2}^{n}\in{\mathscr{P}}(\Omega)$ and let $\tau>0$ and $b_{1},b_{2}>0$. Then the set of minimizers of the problem

[TABLE]

is non-empty and any solution $(\rho_{1}^{*},\rho_{2}^{*})\in{\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega)$ is the characteristic function of a partition of $\Omega$.

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-$\star$ compactness of ${\mathscr{P}}(\Omega)\times{\mathscr{P}}(\Omega)$. Let us remark that the constraint $\rho_{1}+\rho_{2}=1$ a.e. is closed under weak convergence, since $\int_{\Omega}\rho_{1}\,{\rm d}x+\int_{\Omega}\rho_{2}\,{\rm d}x=\mathscr{L}^{d}(\Omega).$

To show that an arbitrary solution $(\rho_{1}^{*},\rho_{2}^{*})$ is the characteristic functions of a partition of $\Omega$, let us rewrite equivalently the minimization problem in terms of transport plans $\pi_{i}\in{\mathscr{P}}(\Omega\times\Omega)$. Recall that $\pi_{i}$ is a plan between $\rho_{i}^{n}$ and $\rho_{i}$, whenever

[TABLE]

for any $\varphi,\psi\in C(\Omega).$ Since we are always working with measures $\rho_{i}^{n},\rho_{i}$ that are absolutely continuous w.r.t. $\mathscr{L}^{d}\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}\Omega$, in the new minimization problem below, as we will see, we can restrict our search to plans that have absolutely continuous marginals w.r.t. $\mathscr{L}^{d}\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}\Omega.$ For a measure $\theta\in{\mathscr{P}}(\Omega\times\Omega)$, we use the notation $\theta^{1}:=(P^{x})_{\#}\theta$ and $\theta^{2}:=(P^{y})_{\#}\theta$ to denote its marginals (here $P^{x},P^{y}:\Omega\times\Omega\to\Omega$ stand for the canonical projections from $\Omega\times\Omega$ onto $\Omega$).

Thus, we aim to solve

[TABLE]

Here we denote

[TABLE]

We define moreover

[TABLE]

and

[TABLE]

where we have extended the second marginals of $\pi_{i}$ by 0 outside of $\Omega$.

The minimization is carried out over a weakly compact set and $\mathcal{S}$ is weakly lower semicontinuous and bounded below, thus minimizers exist. If $\bm{\pi}^{*}=(\pi_{1},\pi_{2})$ is a minimizer of (12) then we can construct a minimizer $\bm{\rho}^{*}=(\rho_{1}^{*},\rho_{2}^{*})$ of the original problem by taking $\rho_{i}^{*}=(P^{y})_{\#}\pi^{*}_{i}$.

Now we consider the properties of minimizers. Clearly, $\mathcal{F}$ is Gâteaux differentiable at $\bm{\pi}^{*}$ in the sense that there exists $\delta\mathcal{F}(\bm{\pi}^{*})\in C(\Omega\times\Omega)$ such that

[TABLE]

where $\bm{\pi}+t\bm{\theta}$ is any admissible perturbation of $\bm{\pi}.$ Similarly, as the other terms in the definition of $\mathcal{S}$ are linear in $\pi$, these are in the same way differentiable, therefore $\mathcal{S}$ is Gâteaux differentiable in this sense.

From Lemma 2.2 it follows that $\mathcal{S}$ is concave along line segments $\bm{\pi}+t\bm{\theta}$ ($t\in(-1,1)$), where $\bm{\pi}$ is a feasible point and $\bm{\theta}=(\theta_{1},\theta_{2})$ is a feasible direction at $\bm{\pi}$ i.e.

[TABLE]

and for some $\delta>0$

[TABLE]

in the sense of signed measures, for all $t\in[0,\delta)$. Furthermore, it follows from Lemma 2.2 that $\mathcal{S}$ is strictly concave on line segments $\bm{\pi}+t\bm{\theta}$, if for some $i$ the marginal $\theta_{i}^{2}(y)$ is not [math] for almost every $y$. Therefore, if $\bm{\pi}^{*}=(\pi_{1}^{*},\pi_{2}^{*})$ is a minimizer and $\bm{\theta}$ is a feasible direction at $\bm{\pi^{*}}$ with at least one non-trivial marginal then

[TABLE]

where $\delta\mathcal{S}(\bm{\pi}^{*})$ stands for the first variation of $\mathcal{S}$ at $\bm{\pi}^{*}$ defined as is (13).

Let $\bm{\pi}$ be a feasible solution and let $\rho_{i}(y)=\pi_{i}^{2}(y)$. Now suppose that there exists $0<\alpha<1$ such that the set $\Omega_{\alpha}=\{y\in\Omega:\rho_{1}(y),\rho_{2}(y)\in(\alpha,1-\alpha)\}$ has positive measure. Partition $\Omega_{\alpha}$ into two sets $E_{1},E_{2}$ of equal measure. Then there exist measure preserving maps $T:E_{1}\to E_{2}$ and $S:E_{2}\to E_{1}$ such that $S\circ T$ and $T\circ S$ are the identity almost everywhere on their respective domains (for example one may choose $T=\nabla\psi$ to be the optimal transport map between the densities $\bm{1}_{E_{1}}$ and $\bm{1}_{E_{2}}$ and $S=\nabla\psi^{*}$). Now we construct a feasible direction $\bm{\theta}$ at $\bm{\pi}$ as follows. Let $r(y)=\frac{\rho_{1}(T(y))}{\rho_{2}(y)}$ for $y\in E_{1}$ and we define the signed measures $\theta_{1},\theta_{2}\in{\mathscr{M}}(\Omega\times\Omega)$ as

[TABLE]

i.e.

[TABLE]

and

[TABLE]

i.e.

[TABLE]

for any $\varphi\in C(\Omega\times\Omega).$

Since $\theta_{1}^{2}(y)=\rho_{1}(T(y))>\alpha$ a.e. on $E_{1}$ we see that $\bm{\theta}$ has a nontrivial marginal.

Let us now check that $\bm{\theta}$ is feasible, i.e. it satisfies (14) and (15). If $\beta<\min\{1-\alpha,\alpha\}$ then $\bm{\pi}\pm\beta\bm{\theta}$ defines a non-negative measure. Next, we check that $\bm{\theta}$ satisfies $\theta_{i}^{1}(x)=0$ for a.e. $x\in\Omega$ and $i=1,2$ and $\sum_{i=1}^{2}\theta_{i}^{2}(y)=0$ for a.e. $y\in\Omega$. For $\varphi\in C(\Omega)$, we have

[TABLE]

and

[TABLE]

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 $\varphi\equiv 1$) that

[TABLE]

Now, for $\varphi\in C(\Omega),$ we have

[TABLE]

Note that our arguments in fact show that $-\bm{\theta}$ is also a feasible direction at $\bm{\pi}$. $\pm\bm{\theta}$ have nontrivial marginals, and it is not possible to have both $\langle\delta\mathcal{S}(\bm{\pi}^{*}),\bm{\theta}\rangle>0$ and $-\langle\delta\mathcal{S}(\bm{\pi}^{*}),\bm{\theta}\rangle>0$. Therefore, $\bm{\pi}$ cannot be a minimizer. This allows us to conclude that for any minimizer $\bm{\rho}^{*}$ of the original problem each density $\rho_{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 notation $K_{\varepsilon}:=\sigma\sqrt{\frac{2\pi}{\varepsilon}}G_{\varepsilon}.$

Lemma 2.4**.**

Let $(\rho_{1}^{*},\rho_{2}^{*})$ be an optimizer in (11) and let $(\rho_{1},\rho_{2})$ a pair of probability measures such that $\rho_{1}+\rho_{2}=1$ a.e. Then,

• (i)

we have the following optimality condition

[TABLE]

for a suitable pair of Kantorovich potentials $(\varphi_{1},\varphi_{2})$ in the optimal transport of $\rho_{1}^{*}$ onto $\rho_{1}^{n}$ and $\rho_{2}^{*}$ onto $\rho_{2}^{n}$, respectively.

• (ii)

there exists a function $p:\Omega\to\mathbb{R}$ that is Lipschitz continuous on $\operatorname{spt}(\rho^{*}_{i})$, $i=1,2$ and is such that for any probability densities $\rho_{1},\rho_{2}$ with $\rho_{1}+\rho_{2}=1$ a.e. we have

[TABLE]

moreover, we have

[TABLE]

Proof.

Let $(\rho_{1},\rho_{2})$ be a pair of probability measures such that $\rho_{1}+\rho_{2}=1$ a.e. For $\delta\in[0,1]$ let us consider the competitors $(\rho_{1}^{\delta},\rho_{2}^{\delta}):=(\rho_{1}^{*}+\delta(\rho_{1}-\rho_{1}^{*}),\rho_{2}^{*}+\delta(\rho_{2}-\rho_{2}^{*}))$, which by construction satisfy the constraint.

By optimality, we have

[TABLE]

Using the exact same argument as in [27, Lemma 3.1] to develop the Wasserstein part on the one hand and the first variations of ${\rm{HC}}_{\varepsilon}$ and $\Phi$ on the other hand, we find (16).

For (ii) (similarly as in [21, Proposition 4.7]), let us notice first that (16) can be written in the form

[TABLE]

where $(h_{1},h_{2})\in L^{\infty}(\Omega)\times L^{\infty}(\Omega)$ is such that $\rho_{1}^{*}+\delta h_{1}+\rho_{2}^{*}+\delta h_{2}=1$ and $\rho_{i}^{*}+\delta h_{i}\in[0,1]$ for $i=1,2$. We know from Proposition 2.3 that $\rho_{1}^{*},\rho_{2}^{*}$ forms a partition of $\Omega$. Therefore, we must take $h_{1}\leq 0$ on $\operatorname{spt}(\rho_{1}^{*})$, and $h_{1}\geq 0$ on $\operatorname{spt}(\rho_{2}^{*})$ (and vice-versa for $h_{2}$). To preserve the constraint $\rho_{1}^{*}+\delta h_{1}+\rho_{2}^{*}+\delta h_{2}=1$ and the mass of each density, we must also take $h_{1}+h_{2}=0$ a.e. and $\int_{\Omega}h_{1}\,{\rm d}x=\int_{\Omega}h_{2}\,{\rm d}x=0$.

Now if we set $h_{2}=-h_{1}$, we find that

[TABLE]

for any $h_{1}\in L^{\infty}(\Omega)$ with 0 mean such that $h_{1}\leq 0$ a.e. on $\operatorname{spt}(\rho_{1})$. This implies that there exist constants $C_{1},C_{2}\in\mathbb{R}$ such that

[TABLE]

[TABLE]

From here, we can define the pressure variable as

[TABLE]

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 $\tau\downarrow 0$). By construction, $p$ clearly satisfies the inequality in (17), 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 (18). ∎

2.4. Continuous in time solutions for $\varepsilon>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 $\varepsilon$ 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 $\tau\downarrow 0$. When $\varepsilon$ is held fixed, we have strong compactness on the pressure term. However, similarly to the models in [20, 21], we will lack strong compactness on the density variables, which would mean that in the limit when $\tau\downarrow 0$, only a weaker version of the system will be available (see (25)). Let us also note that later on in Section 3, we will need some of the pressure estimates provided below when we take $\varepsilon$ to zero along with $\tau$.

Let $T>0$ be a given time horizon and $N\in\mathbb{N}$ and $\tau>0$ such that $N\tau=T.$ Let $\varepsilon>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

[TABLE]

Let $\rho_{i}^{\tau}:[0,T]\to\mathcal{P}(\Omega)$ be defined as

[TABLE]

We define the corresponding velocities, pressures and momentum variables as

[TABLE]

where $\varphi_{i}^{n+1}$ and $p_{i}^{n+1}$ are defined in Lemma 2.4. Following the same steps as in [20, 27], the analysis boils down to obtain a sufficient amount of uniform (in $\tau$) estimates and compactness for the previously obtained functions, then pass to the limit $\tau\downarrow 0.$

Also, as additional tools we construct the corresponding continuous in time (geodesic) interpolations between the densities and the corresponding velocities and momentum variables, $(\tilde{\rho}_{i}^{\tau},\tilde{v}_{i}^{\tau},\tilde{E}_{i}^{\tau})$. We refer to [20, Section 3.1] (see also [27, 37]) for the precise construction. In particular, as a consequence of this construction, we have

[TABLE]

in the sense of distribution on $(0,T)\times\Omega$. Let us comment on the role of the geodesic interpolations. These interpolations (beside the more standard piecewise constant ones) in the context of $W_{2}$–gradient flows were first used by Santambrogio (see the discussion in [37, Section 8.3]). Their role is the following: for any $\tau$, these interpolations, since pieces of geodesic curves, by definition solve continuity equations. Since the piecewise constant interpolations match the geodesic ones at node points $n\tau$, if both of them converge, they need to converge to the same limit. Therefore, one automatically has a continuity equation as the limit of piecewise constant interpolations. An alternative way (which is more often used in the literature) would be to say that the piecewise constant interpolations solve a continuity equation up to an error term, then one would need to show that this error term is converging to zero as the time discretization parameter tends to zero.

Based on the same techniques as in [27, 20], it is easy to obtain the following estimates.

Lemma 2.5**.**

Let $(\rho_{i}^{\tau},v_{i}^{\tau},E_{i}^{\tau})$ and $(\tilde{\rho}_{i}^{\tau},\tilde{v}_{i}^{\tau},\tilde{E}_{i}^{\tau})$ be the previously constructed piecewise constant and continuous in time interpolations, respectively. Then there exists $C>0$ independent of $\tau>0$ and depending only on ${\rm{HC}}_{\varepsilon}(\rho_{1,0},\rho_{2,0})$ such that

• (i)

$W_{2}(\rho_{i}^{\tau}(t),\rho_{i}^{\tau}(s))\leq C\sqrt{t-s+\tau}$* and $W_{2}(\tilde{\rho}_{i}^{\tau}(t),\tilde{\rho}_{i}^{\tau}(s))\leq C\sqrt{t-s}$ for any $0\leq s\leq t\leq T.$ Moreover, up to passing to subsequences $(\rho_{i}^{\tau})_{\tau>0}$ and $(\tilde{\rho}_{i}^{\tau})_{\tau>0}$ converge (uniformly with respect to $W_{2}$) as $\tau\downarrow 0$ to the same limit.*

• (ii)

$(v_{i}^{\tau})_{\tau>0}$* is uniformly bounded in $L^{2}([0,T];L^{2}_{\rho_{i}^{\tau}}).$*

• (iii)

$(E_{i}^{\tau})_{\tau>0}$* and $(\tilde{E}_{i}^{\tau})_{\tau>0}$ are uniformly bounded in ${\mathscr{M}}^{d}([0,T]\times\Omega)$ and up to passing to subsequences they have the same distributional limits as $\tau\downarrow 0$.*

• (iv)

$\displaystyle\int_{0}^{T}\int_{\mathbb{R}^{d}}|\nabla G_{\varepsilon}\star\rho_{i}^{\tau}|\,{\rm d}x\,{\rm d}t\leq CT\,{\rm{HC}}_{\varepsilon}(\rho_{1,0},\rho_{2,0}).$**

• (v)

$\displaystyle\int_{0}^{T}\int_{\Omega}|\rho_{i}^{\tau}(t,x+\delta d)-\rho_{i}^{\tau}(t,x)|\,{\rm d}x\,{\rm d}t\leq CT(\delta+\sqrt{\varepsilon}){\rm{HC}}_{\varepsilon}(\rho_{1,0},\rho_{2,0})$* for any $d\in\mathbb{R}^{d}$ with $|d|=1$ and any $\delta>0$.*

• (vi)

$\displaystyle\int_{0}^{T}\int_{\Omega}|\nabla p^{\tau}|^{2}\,{\rm d}x\,{\rm d}t\leq C\int_{0}^{T}\int_{\Omega}\frac{1}{\varepsilon^{2}}\left(\left(G_{2\varepsilon}\star\rho_{2}^{\tau}\right)^{2}\rho_{1}^{\tau}+\left(G_{2\varepsilon}\star\rho_{1}^{\tau}\right)^{2}\rho_{1}^{\tau}\right)\,{\rm d}x\,{\rm d}t+C$, and as a consequence, $(\nabla p^{\tau})_{\tau>0}$ is uniformly bounded in $L^{2}([0,T]\times\Omega;\mathbb{R}^{d})$ and $\left(p^{\tau}-\fint_{\Omega}p^{\tau}(\cdot,x)\,{\rm d}x\right)_{\tau>0}$ is uniformly bounded in $L^{2}([0,T]\times\Omega)$ by a constant of the form $TC(1+1/\varepsilon^{2}).$

• (vii)

$(\nabla p^{\tau})_{\tau>0}$* is uniformly bounded in $L^{2}([0,T];(C^{1}(\Omega))^{*}),$ independently of $\tau$ and $\varepsilon$. In particular, $\nabla p^{\tau}$ is uniformly bounded in $\mathscr{D}^{\prime}((0,T)\times\Omega;\mathbb{R}^{d})$ and $\nabla p^{\tau}(t,\cdot)$ defines a uniformly bounded distribution of order one, for a.e. $t\in[0,T]$.*

Proof.

Let us notice that the proofs of point (i), (ii) and (iii) follow the exact same lines of the proofs of [27, Lemma 3.3] and [20, Lemma 3.3], so we omit them.

From Proposition 2.3 we know that that $\rho^{n}_{i}$’s are characteristic functions, therefore the proofs of (iv) and (v) follow the exact same lines as the proof of [22, Lemma 2.4.].

Let us give the details on (vi). From the identities (18) from Lemma 2.4, we have that there exists $C>0$ (independent of $\tau>0$, which might increase from one inequality to the next) such that

[TABLE]

where we have used the uniform bounds on $(\rho_{i}^{\tau})_{\tau>0}$ from the previous points and

[TABLE]

Since $\rho_{1}^{n}+\rho_{2}^{n}=1$ a.e., this previous bound implies that $(\nabla p^{\tau})_{\tau>0}$ is uniformly bounded in $L^{2}([0,T]\times\Omega;\mathbb{R}^{d}).$ As a consequence of Poincaré’s inequality, we have that $\left(p^{\tau}-\fint_{\Omega}p^{\tau}(\cdot,x)\,{\rm d}x\right)_{\tau>0}$ is uniformly bounded in $L^{2}([0,T]\times\Omega)$. Both uniform bounds have the form $TC(1+1/\varepsilon^{2}).$

To show (vii), let $\xi\in C^{1}([0,T]\times\Omega;\mathbb{R}^{d})$. Fix $t\in[n\tau,(n+1)\tau)$. Taking the inner product of both sides in (18) with $\xi(t,\cdot)$ and integrating on $\Omega$ w.r.t. $\rho_{i}^{n+1}$ (we drop the dependence on $t$ in the notation of $\xi$), we obtain

[TABLE]

and interchanging the roles of $\rho_{1}^{n+1}$ and $\rho_{2}^{n+1}$, we get

[TABLE]

First, we have

[TABLE]

where in the second and third equalities we have used the change of variables $x\mapsto x-\sqrt{\varepsilon}z$ and $z\mapsto-z$, respectively and the fundamental theorem of calculus in the last equality. Now, since

[TABLE]

for some universal constants $\alpha,\beta>0$, by the previous chain of equalities, we obtain

[TABLE]

where, in the last inequality we used the monotonicity of ${\rm{HC}}_{\varepsilon}$ along the sequence $\left(\rho_{1}^{n},\rho_{2}^{n}\right)_{n}$.

Furthermore, we have that

[TABLE]

and

[TABLE]

Using the fact that we have piecewise constant interpolations, integrating the last equality in time on $[0,T]$, we get

[TABLE]

Now, adding together (23) and (24), then integrating in time on $[0,T]$ and using (ii), we obtain

[TABLE]

where the constant $C>0$ is independent of $\tau>0$ and $\varepsilon>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 $\tau\to 0$ when $\varepsilon$ is held fixed.

Theorem 2.1**.**

Let $\varepsilon>0$, $\rho_{1,0},\rho_{2,0}\in{\mathscr{P}}(\Omega)$ such that $\rho_{1,0}+\rho_{2,0}=1$ a.e. There exists $\rho_{i}\in{\rm{AC}}^{2}([0,T];{\mathscr{P}}(\Omega))$, $i=1,2$, $\overline{p}\in L^{2}([0,T];H^{1}(\Omega))$ and $\zeta_{1},\zeta_{2}\in L^{2}([0,T]\times\Omega;\mathbb{R}^{d})$ such that $\rho_{1}+\rho_{2}=1$ a.e. in $[0,T]\times\Omega$, $b_{1}\zeta_{1}+b_{2}\zeta_{2}=\nabla\overline{p}$ a.e. in $[0,T]\times\Omega$ and the system

[TABLE]

is satisfied in the sense of distributions on $(0,T)\times\Omega.$

It remains open whether in the previous theorem we have strong convergence $\rho_{i}^{\tau}\to\rho_{i}$ in $L^{2}([0,T]\times\Omega)$ as $\tau\downarrow 0$, and in particular whether one can claim that $\zeta_{i}=b_{i}^{-1}\rho_{i}\nabla\overline{p}.$

Proof of Theorem 2.1.

Using the uniform (in $\tau$) estimates in Lemma 2.5, we have the existence of $\rho_{i}\in{\rm{AC}}^{2}([0,T];{\mathscr{P}}(\Omega))$, $i=1,2$ such that up to passing to a subsequence, both $(\rho_{1}^{\tau},\rho_{2}^{\tau})$ and $(\tilde{\rho}_{1}^{\tau},\tilde{\rho}_{2}^{\tau})$ converge to them uniformly in time, w.r.t. $W_{2}$.

Since $\left(p^{\tau}-\fint_{\Omega}p^{\tau}(\cdot,x)\,{\rm d}x\right)_{\tau>0}$ is uniformly bounded in $L^{2}([0,T];H^{1}(\Omega))$, there exists $\overline{p}\in L^{2}([0,T];H^{1}(\Omega))$ such that up to passing to a subsequence, $p^{\tau}-\fint_{\Omega}p^{\tau}(\cdot,x)\,{\rm d}x\rightharpoonup\overline{p}$, weakly in $L^{2}([0,T]\times\Omega)$ and $\nabla p^{\tau}\rightharpoonup\nabla\overline{p}$, weakly in $L^{2}([0,T]\times\Omega;\mathbb{R}^{d})$ as $\tau\downarrow 0$.

We only need to identify the limits of the momentum variables $(E_{i}^{\tau})_{\tau>0}$. From the weak convergence of the density variables, we have that up to passing to a subsequence, $\rho_{1}^{\tau}\nabla K_{\varepsilon}\star\rho_{2}^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}\rho_{1}\nabla K_{\varepsilon}\star\rho_{2}$, as $\tau\downarrow 0$ and similarly $\rho_{2}^{\tau}\nabla K_{\varepsilon}\star\rho_{1}^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}\rho_{2}\nabla K_{\varepsilon}\star\rho_{1}$, and $\rho_{i}^{\tau}\nabla\Phi_{i}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}\rho_{i}\nabla\Phi_{i}$, as $\tau\downarrow 0$ as vector measures.

Furthermore, since $(\rho_{i}^{\tau}\nabla p^{\tau})_{\tau>0}$ is uniformly bounded in $L^{2}([0,T]\times\Omega;\mathbb{R}^{d})$, there exists $\zeta_{i}\in L^{2}([0,T]\times\Omega;\mathbb{R}^{d})$, $i=1,2$ such that up to passing to a subsequence $b_{i}^{-1}\rho_{i}^{\tau}\nabla p^{\tau}\rightharpoonup\zeta_{i}$ weakly in $L^{2}([0,T]\times\Omega;\mathbb{R}^{d})$, as $\tau\downarrow 0$.

Last, by the previous arguments, we have

[TABLE]

as $\tau\downarrow 0$ in $L^{2}([0,T]\times\Omega;\mathbb{R}^{d}).$ Therefore, the thesis of the theorem follows.

∎

It is 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 $(\rho_{1}^{\tau},\rho_{2}^{\tau})$ given as above with $\varepsilon>0$ fixed, for any $\alpha\in(0,1)$ we have

[TABLE]

Proof.

Let us set $i=1$, the other case will be parallel. We have

[TABLE]

By Jensen’s inequality, the above is smaller than

[TABLE]

We then estimate

[TABLE]

Since $\rho^{\tau}_{1}$ and $\rho^{\tau}_{2}$ take values in $\{0,1\}$, we have $|\rho_{1}^{\tau}(x)-\rho_{1}^{\tau}(x-\sqrt{\varepsilon}z)|\leq\rho^{\tau}_{1}(x)\rho^{\tau}_{2}(x-\sqrt{\varepsilon}z)+\rho^{\tau}_{1}(x-\sqrt{\varepsilon}z)\rho^{\tau}_{2}(x)$. Applying these inequalities, the result follows. ∎

3. Muskat flow with surface tension

In this section, we complete the proof of Theorem 1.1 and show that when $\varepsilon=\varepsilon_{\tau}$ goes to zero along with $\tau$, the time interpolated minimizing movements scheme constructed in (21)-(22) 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 (18) 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]\times\Omega$. 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 [22].

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 $\Gamma$, the Muskat problem is given by the continuity equation

[TABLE]

along with the free boundary problem for the pressure

[TABLE]

where $\kappa$ is the mean curvature of $\Gamma$, oriented to be positive when $\operatorname{spt}(\rho_{2})$ is convex at the point. The weak formulation of the Muskat problem with surface tension is provided in Definition 3.1 below.

Definition 3.1**.**

We say that $(\rho_{i},v_{i},p)$ is a weak solution to the Muskat problem with surface tension, if for a.e. $T>0$, $\rho_{i}\in L^{1}([0,T];BV(\Omega;\{0,1\}))\cap{\rm{AC}}^{2}([0,T];{\mathscr{P}}(\Omega))$, $\rho_{1}\rho_{2}=0$ and $\rho_{1}+\rho_{2}=1$ a.e., $v_{i}\in L^{2}([0,T];L^{2}_{\rho_{i,t}}(\Omega;\mathbb{R}^{d}))$, $p\in L^{2}([0,T];(C^{0,\alpha}(\Omega))^{*})$ and for any $\psi\in C^{\infty}([0,T]\times\overline{\Omega})$ and for any vector field $\xi\in C^{\infty}([0,T]\times\Omega;\mathbb{R}^{d})$ with zero normal component on $\partial\Omega$, we have

[TABLE]

with

[TABLE]

Remark 3.1*.*

Let us remark here that $\frac{D\rho_{1}}{|D\rho_{1}|}$ stands for the $L^{\infty}$ density of $D\rho_{1}$ with respect to the total variation measure $|D\rho_{1}|$ (after using the Radon-Nikodym differentiation). Furthermore, let us notice that the term $\left(\frac{D\rho_{1}}{|D\rho_{1}|}\otimes\frac{D\rho_{1}}{|D\rho_{1}|}:\nabla\xi\right)\left(|D\rho_{1}|+|D\rho_{2}|\right)$ is meaningful in the sense that $f(\nu/|\nu|)\,{\rm d}|\nu|$ defines a matrix valued element of ${\mathscr{M}}(\Omega)$ for any $f:\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ which is continuous and 1-homogeneous (see for instance [3, Proposition 3.15] in the case when $\nu=D\rho$ for some $\rho\in BV(\Omega)$). In particular, if $\rho\in BV(\Omega;\{0,1\})$, then $D\rho/|D\rho|$ stands for the measure theoretic normal to the boundary of the set $\operatorname{spt}(\rho)$ and by $\nu/|\nu|$ we mean the density of $\nu$ with respect to its total variation measure $|\nu|$.

Remark 3.2*.*

Although we restrict our attention in the weak curvature equation (27) to test functions with vanishing normal component on $\partial\Omega$, we do not lose information at the boundary. The weak continuity equation (26) encodes the zero normal boundary condition for the velocities, and (27) still encodes the condition that $\Gamma$ and $\partial\Omega$ meet orthogonally (c.f. the proof of Lemma 3.1).

Now we are ready to state the main result of our paper.

Theorem 3.1**.**

Given initial data $\rho_{1,0},\rho_{2,0}\in BV(\Omega)$ such that $\rho_{1,0}\rho_{2,0}=0$ and $\rho_{1,0}+\rho_{2,0}=1$ a.e. in $\Omega$ and Lipschitz continuous potential functions $\Phi_{1},\Phi_{2}:\Omega\to\mathbb{R}$, there exists $(\rho_{i},v_{i},p)$ with $i=1,2$ such that under the energy convergence assumption (EC), it solves (MP1)-(MP2) in the sense of Definition 3.1.

We postpone 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 (26) and (27).

Lemma 3.1**.**

If smooth solutions of (MP1)and (MP2) exist with $\rho_{i}=\chi_{A_{i}}$ and with a $C^{2}$ hypersurface $\Gamma=\partial A_{1}=\partial A_{2}$, then $\rho_{i}$ satisfies (26) and (27) with the choice of

[TABLE]

where $\mu\in(C([0,T]\times\Omega))^{*}$ is the surface measure of $\Gamma=\cup_{t>0}(\Gamma_{t}\times\{t\})$, i.e. after disintegration

[TABLE]

or using test functions

[TABLE]

Remark 3.3*.*

Note that our notion of solution requires adding the surface measure $\sigma\,{\rm d}\mu$ 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 (28).

Proof of Lemma 3.1.

(26) is a standard weak expression of (MP1) with $b_{i}v_{i}=-\nabla p_{i}-\nabla\Phi_{i}$. Let us write again $\Gamma=\cup_{t>0}(\Gamma_{t}\times\{t\})$. Then we have

[TABLE]

Here for the second equality we used integration by parts for the first integral, using the fact that

[TABLE]

and for the third equality we used the curvature jump condition at the interface, and the fact that $\xi$ has zero normal component on $\partial\Omega$. To conclude, let us recall that $\kappa$ is oriented to be positive when $\operatorname{spt}(\rho_{2})$ is convex at the point. With this orientation, observe that for any $C^{1}$ vector field $\xi$ we have

[TABLE]

where $\nu=D\rho_{1}/|D\rho_{1}|$ is normal to $\Gamma$ toward the support of $\rho_{1}$, $\tilde{n}$ stands for the co-normal vector (orthogonal to $\partial\Gamma_{t}$ and tangential to $\Gamma_{t}$). Note that the lower dimensional term $\int_{\partial\Gamma_{t}\cap\partial\Omega}\xi\cdot\tilde{n}\,{\rm d}{\mathscr{H}}^{d-2}$ must vanish since we know that the co-normal $\tilde{n}$ coincides with the boundary normal $n$. Indeed, we can write

[TABLE]

where the final equality follows from the fact that $\xi$ has zero normal component on $\partial\Omega$.

∎

3.1. Preliminary estimates

We present below a compactness result on the piecewise constant interpolations of the density variables $(\rho_{1}^{\tau},\rho_{2}^{\tau})_{\tau>0}$, when the parameter $\varepsilon$ is vanishing together with $\tau.$

Proposition 3.2**.**

Let $(\rho_{1}^{\tau},\rho_{2}^{\tau})_{\tau>0}$ be the piecewise constant interpolations constructed in Section 2 using the minimizing movements scheme (20). Let moreover $\rho_{1,0},\rho_{2,0}\in{\mathscr{P}}(\Omega)$ be such that $\rho_{i,0}\in BV(\Omega;\{0,1\})$ with $\rho_{1,0}+\rho_{2,0}=1$ a.e. in $\Omega$ and in particular $\mathcal{E}(\rho_{1,0},\rho_{2,0})<+\infty.$

Then, there exists $\rho_{1},\rho_{2}\in L^{1}([0,T];BV(\Omega;\{0,1\}))$ such that $\rho_{1}+\rho_{2}=1$ a.e. in $[0,T]\times\Omega$ and up to passing to a subsequence, $\rho_{i}^{\tau}\to\rho_{i}$ strongly in $L^{1}([0,T]\times\Omega)$ as $\varepsilon_{\tau}:=\max\{\tau,\varepsilon\}\downarrow 0$.

Proof.

First, let us notice that by the uniform quasi-Hölder estimate from Lemma 2.5(i), we have that there exists a subsequence of $(\rho_{i}^{\tau})_{\tau>0}$ (that we do not relabel) and $\rho_{i}\in{\rm{AC}}^{2}([0,T];{\mathscr{P}}(\Omega))$ such that $W_{2}(\rho_{i,t}^{\tau},\rho_{i,t})\to 0$ as $\varepsilon_{\tau}\downarrow 0$, uniformly in $t.$ We shall work with this subsequence from now on.

The rest of the proof relies on a careful adaptation of the Aubin-Lions lemma, developed in [36] and used in similar context for instance in [20, 21].

Let us notice that in order to use the Aubin-Lions argument, we need to show a tightness condition (cf. [36, Definition 1.3, Remark 1.5]) for the time-dependent family $(\rho_{i}^{\tau})_{\tau>0}$. This is typically done by providing a compact set (of the space of probability measures), where ‘most’ of the sequences $(\rho_{i}^{\tau}(t))_{\tau>0}$ lie. Note that the estimate in Lemma 2.5(v) does not provide such a compact set. Indeed, for $\tau>0$, the densities $(\rho_{i}^{\tau}(t))_{\tau>0}$ are not actually BV in space. Inspired by arguments in [13], we will work first with an auxiliary sequence $\overline{\rho}_{i}^{\tau}:[0,T]\times\Omega\to[0,1]$, defined as

[TABLE]

where we have performed a convolution with the heat kernel $G_{\varepsilon}$ only in the spacial variable. It worth noticing that this ‘perturbation’ of $\rho_{i}^{\tau}$ is reminiscent to the one obtained via the so-called flow interchange technique introduced in [26]. We have the following properties for this new sequence.

Claim 1. $(\overline{\rho}_{i}^{\tau})_{\tau>0}$ is uniformly bounded (w.r.t. $\varepsilon$ and $\tau$) in $L^{1}([0,T];BV(\Omega)).$

Proof of Claim 1. The uniform $L^{\infty}$ bounds on $(\rho_{i}^{\tau})_{\tau>0}$ are also preserved for $(\overline{\rho}_{i}^{\tau})_{\tau>0}$. Now, as in the proof of [22, Lemma 2.4] we conclude that

[TABLE]

for a uniform constant $C$, which proves our claim.

Claim 2. There exists a subsequence of $(\overline{\rho}_{i}^{\tau})_{\tau>0}$ (that we do not relabel) and $\overline{\rho}_{i}\in L^{1}([0,T]\times\Omega)$ such that $\overline{\rho}_{i}^{\tau}\to\overline{\rho}_{i}$ strongly in $L^{1}([0,T]\times\Omega)$ as $\varepsilon_{\tau}\downarrow 0.$

Proof of Claim 2. For $t\in(0,T)$ fixed, let us notice that $\overline{\rho}_{i,\tau}^{\tau}$ actually can be seen as the evolution of $\rho_{i,t}^{\tau}$ via the heat flow for time $\varepsilon>0$. It is well-known that $W_{2}$ is contractive along the heat flow, so we have

[TABLE]

Now let us set $g:L^{1}(\Omega)\times L^{1}(\Omega)\to[0,+\infty]$ and $\mathcal{F}:L^{1}(\Omega)\to[0,+\infty]$ defined as

[TABLE]

and

[TABLE]

By construction, $\mathcal{F}$ is convex, l.s.c. in $L^{1}(\Omega)$ and it sublevel sets are compact in $L^{1}(\Omega),$ therefore it defines a normal coercive integrand.

From (29) and from the definition of $\mathcal{F}$, we have

[TABLE]

And similarly from Lemma 2.5(i) and from (30), we have

[TABLE]

therefore the assumptions of [36, Theorem 2] are fulfilled and one can conclude that there exists $\overline{\rho}_{i}\in L^{1}([0,T]\times\Omega)$ and a subsequence of $(\overline{\rho}_{i}^{\tau})_{\tau>0}$ (that we do not relabel) which is converging in measure, and in particular pointwise a.e. to $\overline{\rho}_{i}$ as $\varepsilon_{\tau}\downarrow 0$. The strong convergence in $L^{1}([0,T]\times\Omega)$ follows from Lebesgue’s dominated convergence theorem, since $(\overline{\rho}_{i}^{\tau})_{\tau>0}$ is uniformly bounded. The claim follows.

Claim 3. There exists a subsequence of the original sequence $(\rho_{i}^{\tau})_{\tau>0}$ which is converging to $\overline{\rho}_{i}$ strongly in $L^{1}([0,T]\times\Omega)$ as $\varepsilon_{\tau}\downarrow 0$.

Proof of Claim 3. Passing to the same subsequence in the original sequence $(\rho_{i}^{\tau})_{\tau>0}$ (that we do not relabel) as in the last step of the proof of Claim 2, we have

[TABLE]

for a constant $C>0$ (independent of $\tau$ and $\varepsilon$) and the claim follows by taking $\varepsilon_{\tau}\downarrow 0.$ In the last inequality we have used

[TABLE]

where we relied on the fact that since $\rho^{\tau}_{1}$ and $\rho^{\tau}_{2}$ take values in $\{0,1\}$, we have $|\rho_{1}^{\tau}(x)-\rho_{1}^{\tau}(x-\sqrt{\varepsilon}z)|\leq\rho^{\tau}_{1}(x)\rho^{\tau}_{2}(x-\sqrt{\varepsilon}z)+\rho^{\tau}_{1}(x-\sqrt{\varepsilon}z)\rho^{\tau}_{2}(x)$.

To conclude, let us notice that since $(\rho_{i}^{\tau})_{\tau>0}$ converges uniformly w.r.t. $W_{2}$ to $\rho_{i}$, $\overline{\rho}_{i}$ and $\rho_{i}$ have to coincide. Also, for this last subsequence, we can pass to the limit as $\varepsilon_{\tau}\downarrow 0$ in the estimation of Lemma 2.5(v) to obtain that actually $\rho_{i}\in L^{1}([0,T];BV(\Omega;\{0,1\})).$

∎

Let $\bm{\rho}=(\rho_{1},\rho_{2})$ the limit point obtained in Proposition 3.2. Later in this section, we will profit from the assumption

[TABLE]

3.2. Derivation of the weak curvature equation for $\varepsilon$ going to zero together with $\tau$

Let $\xi\in C^{2}([0,T]\times\Omega;\mathbb{R}^{d})$. Fix $t\in[n\tau,(n+1)\tau).$ We consider the piecewise constant interpolations $(\rho_{i}^{\tau},v_{i}^{\tau},p^{\tau})$. Let us take the inner product of the equations (18) with $\xi(t,\cdot)$, multiply with the corresponding $\rho_{i}^{\tau}$ and integrate over $[0,T]\times\Omega$. Adding up the two equations we get

[TABLE]

rearranging the terms yields

[TABLE]

Our aim is now to pass to the limit in this last expression (31) as $\varepsilon_{\tau}\downarrow 0$ to recover (27).

Theorem 3.2**.**

Let $(\rho_{i}^{\tau},v_{i}^{\tau},p^{\tau})_{\tau>0}$ be the piecewise constant interpolations constructed in (21)-(22). Then there exists $\rho_{i}\in L^{1}([0,T];BV(\Omega;\{0,1\}))\cap{\rm{AC}}^{2}([0,T];{\mathscr{P}}(\Omega))$, $v_{i}\in L^{2}([0,T];L^{2}_{\rho_{i}}(\Omega;\mathbb{R}^{d}))$ and $p\in{L^{2}([0,T];(C^{0,\alpha}(\Omega))^{*})}$ such that, up to passing to a subsequence that we do not relabel, we have the following

• (i)

$\rho_{i}^{\tau}\to\rho_{i}$, as $\varepsilon_{\tau}\downarrow 0$, strongly in $L^{1}([0,T]\times\Omega)$;

• (ii)

$v_{i}^{\tau}\rho_{i}^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}v_{i}\rho_{i}$, as $\varepsilon_{\tau}\downarrow 0$, weakly-$\star$* in ${\mathscr{M}}^{d}([0,T]\times\Omega)$;*

• (iii)

$\nabla p^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}\nabla p$, as $\varepsilon_{\tau}\downarrow 0$, weakly-$\star$* in $L^{2}([0,T];(C^{1}(\Omega))^{\star});$*

Proof.

(i) is a consequence of Proposition 3.2 via the Aubin-Lions type argument.

(ii) Let us notice that the estimates in Lemma 2.5(i-iii) are independent of $\varepsilon>0$, therefore there exists $E_{i}\in{\mathscr{M}}^{d}([0,T]\times\Omega)$, $i=1,2$ such that $E_{i}^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}E_{i}$ and $\tilde{E}_{i}^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}E_{i}$ as $\varepsilon_{\tau}\downarrow 0$. Moreover, we have that $(\rho_{i},E_{i})$ solves $\partial_{t}+\nabla\cdot(E_{i})=0$ in the sense of distributions and $\rho_{i}\in{\rm{AC}}^{2}([0,T];({\mathscr{P}}(\Omega),W_{2})).$ Therefore, there exists $v_{i}\in L^{2}([0,T];L^{2}_{\rho_{i}}(\Omega;\mathbb{R}^{d}))$ such that $\partial_{t}\rho_{i}+\nabla\cdot(\rho_{i}v_{i})=0$ in the sense of distributions.

(iii) Let us notice first that from Lemma 2.5(vii) we have that the sequence $(\nabla p^{\tau})_{\tau}$ is uniformly bounded in $L^{2}([0,T];(C^{1}(\Omega))^{*})$, independently of $\varepsilon$. Then, the Banach-Bourbaki-Alaouglu theorem yields that it is sequentially pre-compact in that space, so we have that there exists $\zeta\in L^{2}([0,T];(C^{1}(\Omega))^{*})$ such that up to passing to a subsequence $\nabla p^{\tau}\stackrel{{\scriptstyle\star}}{{\rightharpoonup}}\zeta$, as $\varepsilon_{\tau}\downarrow 0$.

Thus, it only remains to show that $\zeta$ is a gradient. Let us notice that by construction of $p^{\tau}$,

[TABLE]

for any incompressible smooth field $\xi$. Therefore, we also have that $\langle\zeta,\xi\rangle=0$ for any incompressible smooth field $\xi$. So we would be done, if one would have a Helmholtz decomposition in this corresponding space. This result is well known (see for instance [39, Lemma 2.2.1]), if $\zeta(t,\cdot)\in W^{-1,q}(\Omega)^{d}$, for some $q\in(1,+\infty)$. However, our limit object has slightly worse regularity.

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 [15, Theorem 5.23-Theorem 5.24])

Claim. Let $\varphi\in C^{0,\alpha}(\Omega)$. Then the problem $-\Delta u=\varphi$ (with homogeneous zero Neumann boundary condition, if $\int_{\Omega}\varphi\,{\rm d}x=0$) has a unique (modulo constants) solution $u\in C^{2,\alpha}(\Omega)$ such that $\|u\|_{C^{2,\alpha}}\leq\|\varphi\|_{C^{0,\alpha}}.$

Now, let $\varphi\in C^{0,\alpha}(\Omega)$ and $u\in C^{2,\alpha}(\Omega)$ as in the claim. Let $t\in[0,T]$. Then we have

[TABLE]

where in the first inequality we have used the last uniform estimate from the proof of Lemma 2.5(vii). This implies that $(p^{\tau}(t,\cdot))_{\tau}$ is uniformly bounded in $(C^{0,\alpha}(\Omega))^{*}$.

To conclude the thesis of the lemma, we observe that (by possibly subtracting the mean) $(p^{\tau})_{\tau}$ is also converging weakly-$\star$ to some $p$, therefore we will have that $\zeta=\nabla p$ in $\mathscr{D}^{\prime}((0,T)\times\Omega).$

∎

To complete the proof of of Theorem 1.1, it remains to show that limit of equation (31) converges to the weak formulation of the Muskat problem. This result is proven in Proposition A.1, which in turn uses the results of Lemmas A.2 and A.3. These are direct consequences of the corresponding results from [22]. However, for completeness and to facilitate the reading, we collected them in the Appendix A below. These results allow to simply conclude this section with the proof of Theorem 3.1.

Proof of Theorem 3.1.

This proof is a direct consequence of Theorem 3.2 and Proposition A.1. Indeed, these two result allow us to pass to the limit in the equation (31) to obtain (27). ∎

4. Numerics and equilibrium shapes

In this section we present several examples of the performance of our numerical scheme and a discussion of equilibrium shapes.

4.1. Numerical implementation

The Muskat problem evolution can be simulated by discretizing the minimizing movements scheme (9) onto a regular grid. At first glance, numerically solving the discretized variational problem is not so simple. Indeed, Problem (9) is not convex with respect to $\rho$ and the Wasserstein distance is challenging to work with numerically. However, as noted in the introduction, the scheme can be substantially simplified by applying the heat content linearization trick used in [13]. To that end, note that the convexity of the heat content gives us the upper bound

[TABLE]

where the second term is the linearization of the heat content about the previous iterate $\bm{\rho}^{n}$. Thus, if we replace (9) by the linearized scheme,

[TABLE]

we obtain a convex variational problem, and inequality (32) ensures that the scheme still posses the energy dissipation property

[TABLE]

A nearly identical argument to the proof of Proposition 2.3 shows that the set of minimizers of (33) scheme always contains a configuration where $\rho_{1}$ and $\rho_{2}$ are characteristic functions.

To solve problem (33) we introduce the pressure as a Lagrange multiplier for the incompressibility constraint, and instead work with the corresponding dual problem. Up to a constant, the dual problem has the form

[TABLE]

where

[TABLE]

and $c_{1}$, $c_{2}$ denote the quadratic $c$-transform

[TABLE]

(note $c$-transforms play an essential role in optimal transport see for instance [37]). The dual problem is concave with respect to $p$, and can be solved using the recently introduced back-and-forth method [18], which efficiently solves optimal transport problems in dual form.

Due to the two phase nature of the problem, the optimal densities $\rho_{i}^{n+1}$ are not a simple function of the optimal pressure $p^{n+1}$ (this is in contrast to one-phase incompressible fluid flow where the occupied region is the support of the pressure). On the other hand, once we have solved for the optimal pressure in (34), we can recover the velocities $v_{i}$ for each phase (c.f. equation (31)). Thus, in principle, one can recover the densities $\rho^{n+1}_{i}$ from $v_{i}$ and $\rho_{i}^{n}$ by solving the continuity equation for time $\tau$. However, solving the continuity equation accurately is challenging due to the discontinuity of the densities at the phase boundary. Luckily, since we know that the densities remain as characteristic functions, we can instead compute $\rho_{i}^{n+1}$ using the level set method [33]. If we let $\varphi$ be the signed distance function to the interface between $\rho_{1}^{n}$ and $\rho_{2}^{n}$, then by solving the transport equation

[TABLE]

for time $\tau$, we can recover $\rho_{i}^{n+1}$ through the sign of $\varphi$. The advantage of this approach is that the transport equation with Lipschitz initial data can be solved much more accurately than the continuity equation with discontinuous initial data.

4.2. Numerical Experiments

We demonstrate the performance of the numerical scheme on 3 different examples in 2 dimensions. In each experiment, we take our computational domain to be the unit square $[0,1]^{2}$ and set the surface tension constant to be $\sigma=0.15$.

In the first two examples, shown in Figures 1 and 2, we and choose potentials $\Phi_{i}(x,y)=-w_{i}y$ where $w_{1}=5$ and $w_{2}=1$. In Figure 1, the starting configuration for phase 1 is a small square and in Figure 2, the starting configuration for phase 1 is a large square. In both cases, the square becomes round and falls to the bottom of the computational domain. However, due to the difference in mass between examples 1 and 2, the equilibrium configurations are different. In Figure 1, the equilibrium configuration is a half disc sitting at the bottom of the domain, while in Figure 2 the equilibrium configuration is a flat strip.

In the last example, shown in Figure 3, we choose a different potential that leads to a topological change. We set

[TABLE]

and $\Phi_{2}(x,y)=0$. The potential encourages phase 1 to migrate to the top and the bottom of the computational domain, with a stronger force attracting the drop to the bottom. Because the potential pulls the drop in opposite directions, ultimately the initial drop is ripped apart into two separate droplets. Thanks to the scheme’s implicit representation of the interface $\Gamma$, there is no difficulty in simulating topological changes.

4.3. Discussions on the structure of equilibrium shapes

In this subsection we discuss global equilibrium of the energy with gravity potentials

[TABLE]

where $\Phi_{i}(x)=c_{i}x_{d}$, with $0<c_{1}<c_{2}$. The order of the constants denote that $\rho_{1}$ is the lighter fluid, where the vector $-e_{d}$ denotes the direction of gravity. The coordinate here is $x=(x^{\prime},x_{d})$ and for simplicity we consider a cylindrical domain,

[TABLE]

Here the convolution is taken with the extension of the density functions as zero outside of the domain, with the density constraint $\rho_{1}+\rho_{2}=1$ in $\Omega$ and $\int_{\mathbb{R}^{d}}\,{\rm d}\rho_{i}=M_{i}$. Since the densities are extended by zero outside $\Omega$, this will produce a Neumann boundary condition for the interface $\Gamma$. In particular, we expect that $\Gamma$ intersects $\partial\Omega$ orthogonally.

Let us mention that away from the global equilibrium, there are diverse possibilities of stationary states for $\tilde{\mathcal{E}}_{\varepsilon}$ even with zero potentials. For instance any choice of characteristic functions $\rho_{1}$ and $\rho_{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 $\varepsilon\to 0$, the $\Gamma$-convergence properties indicate that the global equilibrium of the $\varepsilon$-energy converges to the limiting density pair $(\rho_{1},\rho_{2})=(\chi_{A},\chi_{A^{c}})$, which is the global minimizer of the limit energy

[TABLE]

under the volume constraint $\mathscr{L}^{d}(A)=M$. Away from the domain boundary, the classical minimal surface theory yields the $C^{2,\alpha}$ regularity of $\partial A$ with the Euler-Lagrange equation

[TABLE]

where $\lambda$ is the Lagrange multiplier associated to the volume constraint and $\kappa$ stands for the curvature.

When $\partial 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 [16]). When the droplet boundary touches the cylinder boundary, various shapes of drops are possible (see Figure 4) and the complete description of possibly non-smooth global minimizers appear to be open. In general, when $\Sigma$ has $C^{1,\alpha}$ boundary, it is shown in [41] that $\partial A$ is $C^{1,\alpha}$ up to the boundary and meets $\partial\Sigma$ orthogonally. For further discussion of available results we refer to [25]. An ongoing work on numerical simulations for our flow suggest that the equilibrium states even for the $\varepsilon$-energy can be categorized as the ones on Figure 4. In dimension two, we could observe an additional equilibrium shape, as in Figure 5.

Acknowledgement. We thank Tim Laux for the helpful discussions and references. We thank Hugo Lavenant for his useful remarks, which in particular led us to fix a small issue in the proof of Proposition 2.3. We thank the two anonymous referees for their very careful reading of the manuscript. Their large amount of comments and remarks helped us to improve the manuscript significantly.

I.K. is partially supported by NSF DMS-1566578. M.J. is partially supported by Simons Math + X Investigator - 510776, DARPA FA8750-18-2-0066, and NSF ATD-1737770. A.R.M was partially supported by the Air Force under the grant AFOSR MURI FA9550-18-1-0502.

Appendix A Passing to the limit in the weak curvature equation based on [22]

Below we collected the results from [22] that are used in the proof of Theorem 3.1. For two of the results we present their proofs as well, either because there was a minor difference between the result that we need and the one from [22] or because we have found a shorter proof than the one in [22].

Let us note that the energy convergence assumption in Theorem 3.1 is used only in Lemma A.2. 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 A.1**.**

For a smooth rapidly decaying kernel $J:\mathbb{R}^{d}\to\mathbb{R}$ and a vector valued Radon measure $\nu\in{\mathscr{M}}^{d}(\Omega)$ we define the Radon measure $\sigma_{J}(\nu)\in{\mathscr{M}}(\Omega)$ as

[TABLE]

Proposition A.1**.**

Let $(\rho_{i}^{\tau})_{\tau>0}$ be the piecewise constant interpolations constructed in (21)-(22) and let $\rho_{i}\in L^{1}([0,T];BV(\Omega;\{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

[TABLE]

${\rm{as\ }}\varepsilon_{\tau}\downarrow 0$, for any $\xi\in C^{3}([0,T]\times\Omega;\mathbb{R}^{d})$. Here, we denoted $J(z):=\sigma\sqrt{2\pi}|z\cdot e_{1}|^{2}G(z)$ and we used the decomposition $\nabla\xi^{{\rm{sym}}}(x)=\sum_{k}\zeta_{k}(x)n_{k}\otimes n_{k}$ with $\zeta_{k}\in C^{\infty}(\Omega)$ and $n_{k}\in S^{d-1}$ such that $\{n_{k}\otimes n_{k}\}_{k=1}^{d(d+1)/2}$ is an appropriate basis of the the space of symmetric matrices.

Proof.

Let $\xi\in C^{3}([0,T]\times\Omega;\mathbb{R}^{d})$. Let us show that along a subsequence

[TABLE]

${\rm{as\ }}\varepsilon_{\tau}\downarrow 0$. This result is not straight forward, since both the functional and the densities depend on $\tau$. This difficulty is the key reason that we need the energy convergence assumption (EC).

Arguing exactly as in the proof of [22, Lemma 2.8], in order to show the convergence result (35), it is enough to show its time-independent version, i.e.

[TABLE]

${\rm{as\ }}\varepsilon_{\tau}\downarrow 0$, for $\mathscr{L}^{1}$-a.e. $t\in[0,T].$

A computation similar to the one in the proof of Lemma 2.5(vii) reveals

[TABLE]

where, we were using a second order Taylor expansion (and the smoothness of $\xi$) in the last equality. Let us observe that the very last term in (36) is converging to 0 as $\varepsilon_{\tau}\downarrow 0$ (due to the strong convergence $\rho_{i}^{\tau}\to\rho_{i}$ in $L^{1}(\Omega)$ and the fact that $\rho_{1}\rho_{2}\equiv 0$ a.e., see Proposition 3.2).

Therefore, it remains to prove

[TABLE]

as $\varepsilon_{\tau}\downarrow 0$. We note that

[TABLE]

where $\nabla\xi^{\textrm{sym}}$ denotes the symmetric part of $\nabla\xi$. A basis for the space of $d\times d$ symmetric matrices is given by the $d+\frac{d(d-1)}{2}$ matrices

[TABLE]

where $e_{i}\in\mathbb{R}^{d}$ is the $i^{th}$ standard basis vector. All of these basis matrices have the form $n\otimes n$ for some $n\in S^{d-1}$. Thus, we may write

[TABLE]

where $\{n_{k}\}$ is any indexing of the above matrices and $\zeta_{k}(x)n_{k}\otimes n_{k}=P_{k}\nabla\xi^{\textrm{sym}}(x)$ where $P_{k}$ is the appropriate projection matrix.

Therefore we have

[TABLE]

Since the functions $z\mapsto|z\cdot n_{k}|^{2}G(z)$ are rotation invariant, when integrating on $\mathbb{R}^{d}$, it is enough to consider only $n_{k}=e_{1}$, where $e_{1}$ is the first element of the standard basis on $\mathbb{R}^{d}$. Now, defining $J(z):=\sigma\sqrt{2\pi}|z\cdot e_{1}|^{2}G(z)$, we deduce the claim from Lemma A.2 and Lemma A.3.

Now, let us remark that a computation completely parallel to [22, Proof of Lemma 3.6.] reveals furthermore that

[TABLE]

as $\varepsilon_{\tau}\to 0$, and so, in this case we would have

[TABLE]

This would in fact complete the claim, as we can compute

[TABLE]

where we have used $\operatorname{Tr}(n_{k}\otimes n_{k})=1$ and the fact that the antisymmetric parts of $\nabla\xi$ are annihilated by the above operations. ∎

The conclusion in the previous proposition is made by the following two lemmas.

Lemma A.2**.**

If $\bm{\rho}^{\tau}\to\bm{\rho}$ in $L^{1}(\Omega)\times L^{1}(\Omega)$ and ${\rm{HC}}_{\varepsilon_{\tau}}(\bm{\rho}^{\tau})\to\mathcal{E}(\bm{\rho})$ as $\varepsilon_{\tau}\downarrow 0$, and $J$ is a nonnegative kernel with rapid decay (i.e. $J(z)\leq|P(z)|G(z)$ for some polynomial $P$) then

[TABLE]

Proof.

The proof of this result is the same as the one of [22, Lemma 3.7]. ∎

Thanks to Lemma A.2 we just need the following pointwise convergence result. This applies to the multiphase case and is stronger than what is needed here. Similarly to [22, Lemma 2.8, Lemma 3.6] we can formulate the following result.

Lemma A.3**.**

Suppose that $\rho_{1},\rho_{2}\in\rm{BV}(\Omega;\{0,1\})$ such that $\rho_{1}(x)\rho_{2}(x)=0$ for a.e. $x\in\Omega$. Then for any smooth function $\zeta:\Omega\to\mathbb{R}$ and any even nonnegative kernel $J:\mathbb{R}^{d}\to\mathbb{R}$ with rapid decay we have

[TABLE]

where we define the Radon measure $\sigma_{J}(D\rho_{i})\in{\mathscr{M}}(\Omega)$ as in Definition A.1 and we have used the notation $J_{\varepsilon}(z):=J(z/\varepsilon).$

The proof supplied below is different from the one by Laux-Otto in [22, Lemma 2.8, Lemma 3.6]. Instead of disintegrating on $\mathbb{R}^{d}$ and using one-dimensional arguments we obtain upper and lower bounds using mollifiers.

Proof.

We begin by showing

[TABLE]

which amounts to showing that

[TABLE]

Expanding out the convolution we have

[TABLE]

Changing variables $x\mapsto x-\varepsilon z$ and then $z\mapsto-z$ in the second term of the integral we get

[TABLE]

Thus, the dominated convergence theorem with the fact that $\rho_{1}\rho_{2}=0$ a.e. yield that the quantity vanishes as $\varepsilon\to 0$.

Now we can restrict our attention to the limit

[TABLE]

Since $\rho_{i}\in\rm{BV}(\Omega;\{0,1\})$ and $\rho_{1}(x)\rho_{2}(x)=0$ a.e.,

[TABLE]

which follows from directly by evaluating both sides. Thus, it suffices to prove

[TABLE]

for any $\chi\in\rm{BV}(\Omega;\{0,1\})$.

For $\delta>0$ let $\eta_{\delta}$ be a smooth approximation to the identity, and set $\chi_{\delta}=\eta_{\delta}\star\chi$, such that $\chi_{\delta}\to\chi$, as $\delta\to 0$ in the sense of strict convergence of BV functions (i.e. $\chi_{\delta}\to\chi$ in $L^{1}$ and $\int_{\Omega}|D\chi_{\delta}|\to\int_{\Omega}|D\chi|$ as $\delta\to 0$; cf. [3, Definition 3.14]).

Then, we also have

[TABLE]

Without loss of generality, one may suppose that $\zeta\geq 0$. By Jensen’s inequality, the above is

[TABLE]

Taking $\varepsilon\to 0$ we get the desired upper bound for the limit.

Conversely, for $\delta>0$ fixed we have from Jensen’s inequality

[TABLE]

Taking $\varepsilon\to 0$ we get

[TABLE]

Finally,

[TABLE]

The result follows.

∎