1 Introduction

# A gradient flow approach to a thin film approximation of the Muskat problem

## Abstract.

A fully coupled system of two second-order parabolic degenerate equations arising as a thin film approximation to the Muskat problem is interpreted as a gradient flow for the -Wasserstein distance in the space of probability measures with finite second moment. A variational scheme is then set up and is the starting point of the construction of weak solutions. The availability of two Liapunov functionals turns out to be a central tool to obtain the needed regularity to identify the Euler-Lagrange equation in the variational scheme.

###### Key words and phrases:
thin film, degenerate parabolic system, gradient flow, Wasserstein distance
###### 2010 Mathematics Subject Classification:
35K65, 35K40, 47J30, 35Q35
Partially supported by the french-german PROCOPE program 20190SE

## 1. Introduction

The Muskat model is a free boundary problem describing the motion of two immiscible fluids with different densities and viscosities in a porous medium (such as intrusion of water into oil). Assuming that the thickness of the two fluid layers is small, a thin film approximation to the Muskat problem has been recently derived in [10] for the space and time evolution of the thickness and of the two fluids ( being then the total height of the layer) and reads

 {∂tf=(1+R)∂x(f∂xf)+R∂x(f∂xg),∂tg=Rμ∂x(g∂xf)+Rμ∂x(g∂xg),(t,x)∈(0,∞)×R, (1.1a) supplemented with the initial conditions f(0)=f0,g(0)=g0,x∈R. (1.1b)

Here, and are two positive real numbers depending on the densities and the viscosities of the fluids. Since and may vanish, (1.1a) is a strongly coupled degenerate parabolic system with a full diffusion matrix due to the terms and . There is however an underlying structure which results in the availability of an energy functional

 E(f,g):=12∫R[f2+R(f+g)2]dx, (1.2)

which decreases along the flow. More precisely, a formal computation reveals that

 ddtE(f,g)=−∫R[f ((1+R)∂xf+R∂xg)2+RRμ g (∂xf+∂xg)2]dx. (1.3)

A similar property is actually valid when (1.1a) is set on a bounded interval with homogeneous Neumann boundary conditions: in that setting, the stationary solutions are constants and the principle of linearized stability is used in [10] to construct global classical solutions which stay in a small neighbourhood of positive constant stationary states. Local existence and uniqueness of classical solutions (with positive components) are also established in [10] by using the general theory for nonlinear parabolic systems developed in [2]. Weak solutions have been subsequently constructed in [9] by a compactness method: the first step is to study a regularized system in which the cross-diffusion terms are “weakened” and to show that it has global strong solutions, the proof combining the theory from [2] for the local well-posedness and suitable estimates for the global existence. Some of these estimates turn out to be independent of the regularisation parameter and provide sufficient information to pass to the limit as the regularisation parameter goes to zero and obtain a weak solution to (1.1a) in a second step. A key argument in the analysis of [9] was to notice that there is another Liapunov functional for (1.1a) given by

 H(f,g):=∫R[flnf+RRμ glng]dx, (1.4)

which evolves along the flow as follows:

 ddtH(f,g)=−∫R[|∂xf|2+R |∂xf+∂xg|2]dx.

The basic idea behind the above computation is to notice that an alternative formulation of (1.1a) is

 {∂tf=∂x[f ∂x((1+R)f+Rg)],∂tg=Rμ ∂x[g ∂x(f+g)],(t,x)∈(0,∞)×R,

so that it is rather natural to multiply the -equation by and the -equation by and find nice cancellations after integrating by parts. In this note, we go one step further and observe that a concise formulation of (1.1a) is actually

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂tf=∂x[f ∂x(δEδf(f,g))],RRμ ∂tg=∂x[g ∂x(δEδg(f,g))],(t,x)∈(0,∞)×R, (1.5)

which is strongly reminiscent of the interpretation of second-order parabolic equations as gradient flows with respect to the -Wasserstein distance, see [3, Chapter 11] and [18, Chapter 8]. Indeed, since the pioneering works [12] on the linear Fokker-Planck equation and [15, 16] on the porous medium equation, several equations have been interpreted as gradient flows with respect to some Wasserstein metrics, including doubly degenerate parabolic equations [1], a model for type-II semiconductors [4], the Smoluchowski-Poisson equation [5], some kinetic equations [6, 8], and some fourth-order degenerate parabolic equations [14], to give a few examples, see also [3] for a general approach. As far as we know, the system (1.5) seems to be the first example of a system of parabolic partial differential equations which can be interpreted as a gradient flow for Wasserstein metrics. Let us however mention that the parabolic-parabolic Keller-Segel system arising in the modeling of chemotaxis has a mixed Wasserstein- gradient flow structure [7].

The purpose of this note is then to show that the heuristic argument outlined previously can be made rigorous and to construct weak solutions to (1.1) by this approach. More precisely, let be the convex subset of the Banach space defined by

 K:={h∈L1(R,(1+x2)dx)∩L2(R):h≥0 a.e. and ∫Rh(x)dx=1}, (1.6)

and consider initial data We next denote the set of Borel probability measures on with finite second moment by and the Wasserstein distance on by Recall that, given two Borel probability measures and in

 W22(μ,ν):=infπ∈Π(μ,ν)∫R2|x−y|2dπ(x,y),

where is the set of all probability measures which have marginals and , that is and for all measurable subsets and of Alternatively, is equivalent to

 ∫R2(ϕ(x)+ψ(y))dπ(x,y)=∫Rϕ(x)dμ(x)+∫Rψ(y)dν(y) for all (ϕ,ψ)∈L1(R;R2).

With these notation, our result reads:

###### Theorem 1.1.

Assume that Given and the sequence obtained recursively by setting

 (f0τ,g0τ) := (f0,g0), (1.7) Fnτ(fn+1τ,gn+1τ) := inf(u,v)∈K2Fnτ(u,v), (1.8)

with

 Fnτ(u,v):=12τ (W22(u,fnτ)+RRμ W22(v,gnτ))+E(u,v),(u,v)∈K2, (1.9)

is well-defined. Introducing the interpolation defined by

 fτ(t):=fnτ and gτ(t):=gnτ for t∈[nτ,(n+1)τ) and n≥0, (1.10)

there exist a sequence of positive real numbers, , and functions such that

 (fτk,gτk)⟶(f,g) in L2((0,T)×R;R2) for all T>0. (1.11)

Moreover,

• , ,

• with

and the pair is a weak solution of (1.1) in the sense that

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∫Rf(t) ξdx−∫Rf0 ξdx+∫t0∫Rf(σ)[(1+R)∂xf+R∂xg](σ)∂xξdxdσ=0,∫Rg(t) ξdx−∫Rg0 ξdx+Rμ ∫t0∫Rg(σ)(∂xf+∂xg)(σ)∂xξdxdσ=0, (1.12)

for all and In addition, satisfy the following estimates

 (a) H(f(T),g(T))+∫T0∫R[|∂xf|2+R|∂x(f+g)|2]dxdt≤H(f0,g0), (b) E(f(T),g(T))+12∫T0∫R[f((1+R)∂xf+R∂xg)2+RRμg(∂xf+∂xg)2]dxdt≤E(f0,g0),

for a.e. , and being the functionals defined by (1.2) and (1.4), respectively.

Let us briefly outline the proof of Theorem 1.1: in the next section, we study the variational problem (1.8) and the properties of its minimizers. A key argument here is to note that the availability of the Liapunov functional (1.4) allows us to apply an argument from [14] which guarantees that the minimizers are not only in but also in . This property is crucial in order to derive the Euler-Lagrange equation in Section 2.2. The latter is then used to obtain additional regularity on the minimizers, adapting an argument from [15]. Convergence of the variational approximation is established in Section 3. Finally, three technical results are collected in the Appendix.

As a final comment, let us point out that we have assumed for simplicity that the initial data and are probability measures but that the case of initial data having different masses may be handled in the same way after a suitable rescaling: more precisely, let and denote a solution to (1.1) by . Setting and and recalling that and for all , we realize that solves

 Missing or unrecognized delimiter for \left

with and initial data . The corresponding variational scheme then involves the functional

 12τ (1∥g0∥1 W22(u,F0)+RRμ∥f0∥1 W22(v,G0))+η22 ∥u∥22+R2 ∥∥η u+η−1 v∥∥22,(u,v)∈K2,

to which the analysis performed below (with ) also applies.

## 2. A variational scheme

Given and we introduce the functional

 Fτ(u,v):=12τ (W22(u,f0)+RRμ W22(v,g0))+E(u,v),(u,v)∈K2, (2.1)

and consider the minimization problem

 inf(u,v)∈K2Fτ(u,v). (2.2)

### 2.1. Existence and properties of minimizers

Let us start by proving that, for each , the minimization problem (2.2) has a unique solution in

###### Lemma 2.1.

Given and there exists a unique minimizer of (2.2). Additionally, with

 ∥∂xf∥22+R∥∂x(f+g)∥22≤1τ [H(f0)−H(f)+RRμ (H(g0)−H(g))], (2.3)

where

 H(h):=∫Rhln(h)dxfor h∈L1(R) % such that h≥0 a.e. and hln(h)∈L1(R). (2.4)

Recall that, if , then (see Lemma A.1 below) so that the right-hand side of (2.3) is well-defined.

###### Proof.

The uniqueness of the minimizer follows from the convexity of and and the strict convexity of the energy functional .

We next prove the existence of a minimizer. To this end, pick a minimizing sequence There exists a constant such that

 ∥uk∥2+∥vk∥2 ≤ C,k≥1, (2.5) W2(uk,f0)+W2(vk,g0) ≤ C,k≥1. (2.6)

From (2.5) we obtain at once that there exist and a subsequence of (denoted again by ) such that

 uk⇀fandvk⇀g%in$L2(R).$ (2.7)

Let us first check that Indeed, the nonnegativity of and readily follows from that of and by (2.7) while integrating the inequality with respect to an arbitrary yields

 ∫Ruk(x)x2dx=∫R2x2dπ(x,y) ≤ 2∫R2y2dπ(x,y)+2∫R2|x−y|2dπ(x,y) ≤ 2∫R2f0(y)y2dy+2∫R2|x−y|2dπ(x,y),

which implies by virtue of (2.6) that

 ∫Ruk(x)x2dx≤2∫Rf0(x)x2dx+2W22(uk,f0)≤C,k≥1. (2.8)

Similarly,

 ∫Rvk(x)x2dx≤C,k≥1. (2.9)

Owing to (2.5), (2.8), and (2.9), we deduce from the Dunford-Pettis theorem that and are weakly sequentially compact in We may thus assume (after possibly extracting a further subsequence) that and in whence

 ∫Rf(x)dx=limk→∞∫Ruk(x)dx=1 and ∫Rg(x)dx=limk→∞∫Rvk(x)dx=1.

Finally, combining (2.5), (2.8), and (2.9) with a truncation argument ensure that and both belong to Summarising, we have shown that

The next step is to prove that

 Fτ(f,g)=inf(u,v)∈K2Fτ(u,v).

Indeed, on the one hand, the weak convergence (2.7) implies that

 E(f,g)≤liminfk→∞E(uk,vk).

On the other hand, we recall that the -Wasserstein metric is lower semicontinuous with respect to the narrow convergence of probability measures in each of its arguments, see [3, Proposition 7.1.3], and the weak convergence of in ensures that

 W22(f,f0)≤liminfk→∞W22(uk,f0)% andW22(g,g0)≤liminfk→∞W22(vk,g0).

Consequently,

 Fτ(f,g)≤liminfk→∞Fτ(uk,vk) with (f,g)∈K2,

so that is a minimizer of in .

As a final step, we show that and belong to To this end, we follow the approach developed in [14] and take advantage of the availability of another Liapunov function as already discussed in the Introduction. More precisely, denote the heat semigroup by , that is,

 (Gth)(x):=1√4πt ∫Rexp(−|x−y|24t) h(y)dy,(t,x)∈[0,∞)×R,

for . Since , classical properties of the heat semigroup ensure that for all . Consequently, and we deduce that

 E(f,g) −E(Gtf,Gtg)≤12τ[(W22(Gtf,f0)−W22(f,f0))+RRμ(W22(Gtg,g0)−W22(g,g0))] (2.10)

for all Moreover, for all , we have

 ddtE(Gtf,Gtg)= ∫R[Gtf ∂tGtf+R (Gtf+Gtg) ∂t(Gtf+Gtg)]dx=−∥∂xGtf∥22−R ∥∂xGt(f+g)∥22,

and by integration with respect to time we find that

 1t ∫t0[∥∂xGsf∥22+R∥∂xGs(f+g)∥22]ds≤E(f,g)−E(Gtf,Gtg)tfor all t>0.

Since is non-increasing for we end up with

 ∥∂xGtf∥22+R∥∂xGt(f+g)∥22≤E(f,g)−E(Gtf,Gtg)tfor all t>0. (2.11)

We recall now some properties of the heat flow in connection with the -Wasserstein distance , see [3, 8, 16], these properties being actually collected in [14, Theorem 2.4]. The heat flow is the gradient flow of the entropy functional given by (2.4) for and, for all we have [3, Theorem 11.1.4]

 12ddtW22(Gth,~h)+H(Gth)≤H(~h)for a.e. t≥0. (2.12)

Choosing and in (2.12), we obtain

 12ddt[W22(Gtf,f0)+RRμW22(Gtg,g0)]≤H(f0)−H(Gtf)+RRμ(H(g0)−H(Gtg))

for a.e. Integrating the above inequality with respect to time and using the time monotonicity of and give

 12[W22(Gtf,f0)−W22(f,f0)+RRμ (W22(Gtg,g0)−W22(g,g0))] Missing or unrecognized delimiter for \left ≤t [H(f0)−H(Gtf)+RRμ (H(g0)−H(Gtg))]. (2.13)

Gathering (2.10), (2.11), and (2.13), we find

 ∥∂xGtf∥22+R ∥∂xGt(f+g)∥22≤1τ [H(f0)−H(Gtf)+RRμ (H(g0)−H(Gtg))] (2.14)

for As a direct consequence of (2.14) and the boundedness from below (A.2) of in , and are bounded in and converge to and , respectively, in the sense of distributions as . This implies that both and belongs to and we can pass to the limit as in (2.14) to obtain the desired estimate (2.3) and finish the proof. ∎

### 2.2. The Euler-Lagrange equation

We now identify the Euler-Lagrange equation corresponding to the minimization problem (2.2).

###### Lemma 2.2.

Given and , the minimizer of in satisfies

 ∣∣∣1τ ∫Rξ (f−f0)dx+∫R[((1+R) f ∂xf+R f ∂xg) ∂xξ]dx∣∣∣ ≤ ∥∂2xξ∥∞2τ W22(f,f0), (2.15) ∣∣∣1τ ∫Rξ (g−g0)dx+Rμ ∫R[(g ∂xf+g ∂xg) ∂xξ]dx∣∣∣ ≤ ∥∂2xξ∥∞2τ W22(g,g0), (2.16)

for .

###### Proof.

To derive (2.15)-(2.16) we follow the general strategy outlined in [18, Chapter 8]. According to Brenier’s theorem [18, Theorem 2.12], there are two convex functions and which are uniquely determined up to an additive constant such that

 W22(f,f0)=∫R|x−∂xφ(x)|2 f0(x)dx=infT#f0=f∫R|x−T(x)|2 f0(x)dx, (2.17a) where the infimum is taken over all measurable functions T:R→R pushing f0 forward to f (f=T#f0), i.e. satisfying ∫Bf(x)dx=∫T−1(B)f0(x)dxfor all Borel sets B% of R, and W22(g,g0)=∫R|x−∂xψ(x)|2 g0(x)dx=infS#g0=g∫R|x−S(x)|2 g0(x)dx. (2.17b)

We pick now two test functions and in and define

 fε:=((id+εξ)∘∂xφ)#f0=(id+εξ)#f and gε:=((id+εη)∘∂xψ)#g0=(id+εη)#g (2.18)

for each , where is the identity function on To ease notation we set

 Tε:=id+εξ and Sε:=id+εη, (2.19)

and observe that there is small enough (depending on both and ) such that, for and are diffeomorphisms in . Then, by (2.18), we find the identities

 fε=f∘T−1ε∂xTε∘T−1ε and gε=g∘S−1ε∂xSε∘S−1ε,ε∈(0,ε0]. (2.20)

Observing that and

 ∥fε∥22=∫R|f(x)|2∂xTε(x)dxand∥gε∥22=∫R|g(x)|2∂xSε(x)dx, (2.21)

we clearly have for all and thus . Consequently,

 0≤12τ [W22(fε,f0)−W22(f,f0)+RRμ (W22(gε,g0)−W22(g,g0))]+E(fε,gε)−E(f,g). (2.22)

Concerning the energy , it follows from (2.21) that

 2(E(fε,gε)−E(f,g))=(1+R) Iε1+R Iε2+2R Iε3, (2.23)

with

 Iε1:=∫R(1∂xTε(x)−1) |f(x)|2dx,Iε2:=∫R(1∂xSε(x)−1) |g(x)|2dx, Iε3:=∫R(fε gε−f g)(x)dx.

We now consider the three integrals in the right-hand side of the relation (2.23) separately: since

 Iε1=−ε∫R∂xξ1+ε∂xξ f2dxandIε2=−ε∫R∂xη1+ε∂xη g2dx,

it readily follows from Lebesgue’s dominated convergence theorem that

 limε→0Iε1ε=−∫R∂xξ f2dxandlimε→0Iε2ε=−∫R∂xη g2dx. (2.24)

We next turn to the term involving both and and split it in two terms with

 Iε31:=∫R(fε+f) (gε−g)dxandIε32:=∫R(gε+g) (fε−f)dx.

By (2.20),

 Iε31= ∫R(f∘T−1ε∂xTε∘T−1ε+f)(g∘S−1ε∂xSε∘S−1ε−g)dx =