Philippe Laurençot
email: laurenco@math.univ-toulouse.fr

Bogdan-Vasile Matioc
email: bogdan-vasile.matioc@univie.ac.at

A THIN FILM APPROXIMATION OF THE MUSKAT PROBLEM WITH GRAVITY AND CAPILLARY FORCES

Keywords: Mathematics Subject Classification. 35K65, 35K41, 47J30, 35Q35 thin film, degenerate parabolic system, gradient flow, Wasserstein distance

Existence of nonnegative weak solutions is shown for a thin film approximation of the Muskat problem with gravity and capillary forces taken into account. The model describes the space-time evolution of the heights of the two fluid layers and is a fully coupled system of two fourth order degenerate parabolic equations. The existence proof relies on the fact that this system can be viewed as a gradient flow for the 2-Wasserstein distance in the space of probability measures with finite second moment.

Introduction and main result

The Muskat problem 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). It gives the space and time evolution of the heights f ≥ 0 and g ≥ 0 of the two fluid layers, with the first layer, of height f , located on a impermeable horizontal bottom and the second one, of height g, on top of it, as well as that of the pressure fields inside the fluids. As it involves four unknowns and two free boundaries, one separating the lower and the upper fluid and one separating the upper fluid and the air, it is a complex problem. Recently, using a lubrication approximation, see, e.g., [START_REF] Leal | Advanced transport phenomena. Fluid mechanics and convective transport processes[END_REF]Chapter 5.B], a thin film approximation to the Muskat problem has been derived in [START_REF] Escher | Modelling and analysis of the Muskat problem for thin fluid layers[END_REF] which retains only the heights f and g of the two fluid layers as unknowns and reads

∂ t f =div [f (-A∇∆f -B∇∆g + a∇f + b∇g)] , ∂ t g =div [g (-∇∆f -∇∆g + c∇f + c∇g)] , (t, x) ∈ (0, ∞) × R 2 , (1.1a)
together with initial conditions (f, g)(0) = (f 0 , g 0 ),

x ∈ R 2 .

(1.1b)

The parameters (A, B, a, b, c) involved in (1.1a) depend on the densities, viscosities, and surface tensions of the fluids and are assumed to satisfy (a, b, c) ∈ [0, ∞) 3 , cB = b, and A > B > 0.

(1.2)

Let us recall that the second order terms in (1.1a) account for gravity forces while the fourth order terms result from capillary forces in the original Muskat problem. The problem (1.1a) is a fourth order degenerate parabolic system with a full diffusion matrix. Its parabolicity has been exploited in [START_REF] Escher | Existence and stability of solutions for a strongly coupled system modelling thin fluid films[END_REF] to show the local existence and uniqueness of positive strong solutions to (1.1a) in a bounded interval (0, L) with no slip boundary conditions. Global existence for initial data close to a positive flat steady state is also proved when π 2 (A -B)/L 2 + (a -b) > 0 as well as the stability of this steady state. The existence and stability of flat and non-flat stationary solutions are also discussed according to the values of the parameters. Nonnegative global solutions have also been constructed either when only gravity is taken into account (A = B = 0) and (1.1) is considered in a bounded interval with homogeneous Neumann boundary conditions [START_REF] Escher | Existence and stability of weak solutions for a degenerate parabolic system modelling two-phase flows in porous media[END_REF] or in R [START_REF] Ph | A gradient flow approach to a thin film approximation of the Muskat problem[END_REF], or when gravity forces are discarded (a = b = c = 0) and (1.1) is considered in a bounded interval with no slip boundary conditions [START_REF] Matioc | Nonnegative global weak solutions for a degenerate parabolic system modeling thin films driven by capillarity[END_REF]. An important tool in the above mentioned works is the availability of two Liapunov functionals for (1.1), one being the functional E defined by

E(f, g) := 1 2 R 2 (A -B)|∇f | 2 + B|∇(f + g)| 2 + (a -b)f 2 + b(f + g) 2 dx (1.3)
for (f, g) ∈ H 1 (R 2 ; R 2 ) which decreases along the trajectories of (1.1) according to

d dt E(f, g) = - R 2 f |∇∆(Af + Bg) -∇(af + bg)| 2 + Bg |∇∆(f + g) -b∇(f + g)| 2 dx .
It turns out that there is an underlying structure in (1.1) which allows us to view it as a gradient flow for the energy E defined in (1.3) with respect to the 2-Wasserstein distance in P 2 (R 2 )×P 2 (R 2 ). Recall that P 2 (R 2 ) is the set of Borel probability measures on R 2 with finite second moment and that, given two Borel probability measures µ and ν in P 2 (R 2 ), the 2-Wasserstein distance W 2 (µ, ν) on P 2 (R 2 ) is defined by

W 2 2 (µ, ν) := inf π∈Π(µ,ν) R 4
|x -y| 2 dπ(x, y) ,

where Π(µ, ν) is the set of all probability measures π ∈ P(R 4 ) which have marginals µ and ν, that is, π[U × R 2 ] = µ[U ] and π[R 2 × U ] = ν[U ] for all measurable subsets U of R 2 . This gradient flow structure was actually uncovered in [START_REF] Ph | A gradient flow approach to a thin film approximation of the Muskat problem[END_REF] in the simplified case where the capillary forces are neglected (A = B = 0) and shown to provide a convenient setting to construct weak solutions to (1.1). This is thus the approach we shall use in this paper to prove the existence of nonnegative weak solutions to (1.1), showing additionally the convergence of the variational approximation. It is worth noticing at this point that, if g = 0, the system (1.1) reduces to the single equation

∂ t f = div [f (-A∇∆f + a∇f )] ,
which also has a gradient flow structure with respect to the 2-Wasserstein distance as already shown in [START_REF] Otto | Dynamics of labyrinthine pattern formation in magnetic fluids: a mean-field theory[END_REF][START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF] for A = 0 and [START_REF] Matthes | A family of nonlinear fourth order equations of gradient flow type[END_REF] for a = 0. Let us recall that, since the pioneering works [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF][START_REF] Otto | Dynamics of labyrinthine pattern formation in magnetic fluids: a mean-field theory[END_REF][START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF], several parabolic equations have been interpreted as gradient flows for Wasserstein distances, including the Fokker-Planck equation [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], the porous medium equation [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF], second order nonlocal and/or degenerate parabolic equations [START_REF] Agueh | Existence of solutions to degenerate parabolic equations via the Monge-Kantorovich theory[END_REF][START_REF] Ambrosio | A gradient flow approach to an evolution problem arising in superconductivity[END_REF][START_REF] Blanchet | Convergence of the mass-transport steepest descent scheme for the subcritical Patlak-Keller-Segel model[END_REF], some kinetic equations [START_REF] Carlen | Solution of a model Boltzmann equation via steepest descent in the 2-Wasserstein metric[END_REF][START_REF] Carrillo | Contractions in the 2-Wasserstein length space and thermalization of granular media[END_REF], and fourth order degenerate parabolic equations [START_REF] Lisini | Cahn-Hilliard and thin film equations with nonlinear mobility as gradient flows in weighted-Wasserstein metrics[END_REF][START_REF] Matthes | A family of nonlinear fourth order equations of gradient flow type[END_REF]. Besides (1.1) there does not seem to be many systems of parabolic partial differential equations endowed with a similar structure with the exception of the parabolic-parabolic chemotaxis Keller-Segel system which has a mixed Wasserstein-L 2 gradient flow structure [START_REF] Blanchet | The parabolic-parabolic Keller-Segel system with critical diffusion as a gradient flow in R d[END_REF][START_REF] Carrillo | A Wasserstein-l 2 mixed gradient flow approach to the fully parabolic Keller-Segel model[END_REF].

Before stating the main result of this paper, we introduce some notations: let K be the convex subset of P 2 (R 2 ) defined by

K := h ∈ L 1 (R 2 , (1 + |x| 2 )dx) ∩ H 1 (R 2 ) : h ≥ 0 a.e. and h 1 = 1 ,
and set K 2 := K × K. We next recall that the functional H defined by

H(h) := R 2 h ln(h) dx (1.4) is well-defined for h ∈ K, see Lemma A.3.
The main result of this paper is the following:

Theorem 1.1. Assume that (1.
2) is satisfied and let (f 0 , g 0 ) ∈ K 2 . Given τ ∈ (0, 1), we define (f 0 τ , g 0 τ ) := (f 0 , g 0 ) and consider for each n ∈ N the minimization problem inf

(u,v)∈K 2 F n τ (u, v) , (1.5)
where

F n τ (u, v) := 1 2τ W 2 2 (u, f n τ ) + B W 2 2 (v, g n τ ) + E(u, v) , (u, v) ∈ K 2 . For each n ∈ N, the minimization problem (1.5) has a solution (f n+1 τ , g n+1 τ ), which is unique if a ≥ b. Defining the interpolation functions (f τ , g τ ) by (f τ , g τ )(t) := (f n τ , g n τ ) for t ∈ [nτ, (n + 1)τ ) and n ∈ N, there exist a sequence τ k ց 0 and functions (f, g) : [0, ∞) → K 2 such that (f τ k , g τ k )(t) → (f, g)(t) in L 2 (R 2 ; R 2 ) (1.6) for all t ≥ 0. Moreover, (i) (f, g) ∈ L ∞ (0, t; H 1 (R 2 ; R 2 )) ∩ L 2 (0, t; H 2 (R 2 ; R 2 )), (ii) (f, g) ∈ C([0, ∞), L 2 (R 2 ; R 2 )
) and (f, g)(0) = (f 0 , g 0 ), and the pair (f, g) is a weak solution of (1.1) in the sense that

       R 2 (f (t) -f 0 )ξ dx + t 0 R 2 [∆(Af + Bg) div(f ∇ξ) + f ∇(af + bg) • ∇ξ] dx ds = 0 R 2 (g(t) -g 0 )ξ dx + t 0 R 2 [∆(f + g) div(g∇ξ) + cg∇(f + g) • ∇ξ] dx ds = 0 (1.7)
for all t > 0 and ξ ∈ C ∞ 0 (R 2 ). Finally, (f, g) satisfies the following estimates:

H(f (t)) + BH(g(t)) + t 0 D H (f (s), g(s)) ds ≤ H(f 0 ) + BH(g 0 ) , (1.8)
E(f (t), g(t)) + 1 2 t 0 w f (s) 2 2 + B w g (s) 2 2 ds ≤ E(f 0 , g 0 ) (1.9)
for all t > 0, where

D H (f, g) := (A -B) ∆f 2 2 + B ∆(f + g) 2 2 + (a -b) ∇f 2 2 + b ∇(f + g) 2 2
, and w f and w g are two vector fields in L 2 ((0, ∞) × R 2 ; R 2 ) defined as follows: introducing the vector fields (j f , j g ) defined by

j f := -∇(f ∆(Af + Bg)) + ∆(Af + Bg)∇f + f ∇(af + bg) j g := -∇(g∆(f + g)) + ∆(f + g)∇g + cg∇(f + g) in D ′ ((0, ∞) × R 2 ; R 2 ), (1.10)
they actually both belong to L 2 (0, ∞; L 4/3 (R 2 ; R 2 )) and

j f = f w f and j g = √ g w g a.e. in (0, ∞) × R 2 . (1.11)
Notice that Theorem 1.1 provides not only the existence of a weak solution to (1.1) but also the convergence of subsequences of solutions to the variational scheme (1.5) to a solution to (1.1). The proof of Theorem 1.1 thus relies strongly on the study of the minimization problem (1.5) which is performed in Section 2. The availability of the second Liapunov functional H(f ) + BH(g) (which is really a Liapunov functional only if a ≥ b) plays an important role here as it allows us to show by an argument from [START_REF] Matthes | A family of nonlinear fourth order equations of gradient flow type[END_REF] that the minimizers of (1.5) are actually in H 2 (R 2 ). This additional regularity is actually at the heart of the identification of the Euler-Lagrange equation satisfied by the minimizers, see Section 3. The convergence of the scheme and the existence of weak solutions to (1.1) are established in the last section, some care being needed to handle the fourth order terms. Indeed, as indicated in (1.10) and (1.11), due to the degenerate character of the equations, we have to deal with expressions that correspond to the fourth order terms in (1.1a) and are not functions.

Remark 1.2.

(1) In contrast to the dissipation of the energy E, the dissipation D H (f, g) of the functional H(f ) + BH(g) need not be nonnegative and this occurs when a < b, see (1.8). It is yet unclear what information on the dynamics may be retrieved from (1.8) in that case.

(2) For simplicity, we have restricted the analysis to initial data (f 0 , g 0 ) satisfying f 0 1 = g 0 1 = 1. A similar result holds true for arbitrary nonnegative and integrable initial data (f 0 , g 0 ) ∈ H 1 (R 2 ; R 2 ) with finite second moment and may be proved analogously to Theorem 1.1. Introducing (F, G) := (f / f 0 1 , g/ g 0 1 ), this new unknown solves a system with the same structure as (1.1), but with different parameters (depending on f 0 1 and g 0 1 ). (3) Theorem 1.1 is also valid for the one dimensional version of (1.1).

Throughout the paper we use the following notation: x = (x 1 , x 2 ) denotes a generic point of R 2 and the partial derivative with respect to x i is denoted by ∂ i , i = 1, 2. For a sufficiently smooth function h, ∇h := (∂ 1 h, ∂ 2 h) denotes its gradient and D 2 h := (∂ i ∂ j h) 1≤i,j≤2 its Hessian matrix. Finally, Dξ denotes the gradient of a vector field ξ = (ξ 1 , ξ 2 ) ∈ C 1 (R 2 ; R 2 ) and is given by Dξ := (∂ i ξ j ) 1≤i,j≤2 , the divergence of ξ being the trace of Dξ, that is, div(ξ

) = ∂ 1 ξ 1 + ∂ 2 ξ 2 .

The minimization problem

In this section, we show that, given (f 0 , g 0 ) ∈ K 2 and τ > 0, the minimization problem

inf (u,v)∈K 2 F τ (u, v) (2.1)
with the functional F τ defined by

F τ (u, v) := 1 2τ W 2 2 (u, f 0 ) + B W 2 2 (v, g 0 ) + E(u, v) , (u, v) ∈ K 2 , (2.2)
has at least one solution and we study the regularity of the minimizers. Since a -b might be negative, which is the case when the more dense fluid lies on top of the less dense one, cf. [START_REF] Escher | Existence and stability of solutions for a strongly coupled system modelling thin fluid films[END_REF], the first step is to prove that E is bounded from below.

Lemma 2.1. The energy functional E satisfies

E(f, g) ≥ -C 1 + B 2 ∇(f + g) 2 2 + A -B 4 ∇f 2 2 for all (f, g) ∈ K 2 , (2.3)
where C 1 is a positive constant depending only A, B, a, b, and c.

Proof. According to the Gagliardo-Nirenberg-Sobolev inequality there is a positive constant

C 2 > 0 such that h 2 ≤ C 2 h 1/2 1 ∇h 1/2 2
for all h ∈ H 1 (R 2 ).

(2.4) By the definition of K, this inequality yields

h 2 2 ≤ C 2 2 ∇h 2 for all h ∈ K. (2.5)
Owing to (2.5) and Young's inequality, we find, for (f, g)

∈ K 2 , 2E(f, g) =(A -B) ∇f 2 2 + B ∇(f + g) 2 2 + (a -b) f 2 2 + b f + g 2 2 ≥B ∇(f + g) 2 2 + (A -B) ∇f 2 2 -C 2 2 (b -a) + ∇f 2 ≥B ∇(f + g) 2 2 + A -B 2 ∇f 2 2 + A -B 2 ∇f 2 - C 2 2 (b -a) + A -B 2 - C 4 2 (b -a) 2 + 2(A -B) ,
whence (2.3).

We now show the existence of a minimizer to (2.1).

Lemma 2.2. Given (f 0 , g 0 ) ∈ K 2 and τ > 0, there exists a minimizer (f, g) ∈ K 2 of (2.1) which is unique if a ≥ b. Moreover, f and g both belong to H 2 (R 2 ) and

(A -B) ∆f 2 2 + B ∆(f + g) 2 2 + (a -b) ∇f 2 2 + b ∇(f + g) 2 2 ≤ 1 τ [H(f 0 ) -H(f ) + B (H(g 0 ) -H(g))] , (2.6)
the functional H being defined in (1.4).

Proof. Pick a minimizing sequence (u k , v k ) k≥1 of F τ in K 2 . Invoking (2.3) and (2.5) we see that

u k H 1 + v k H 1 ≤ C , k ≥ 1 , (2.7) W 2 (u k , f 0 ) + W 2 (v k , g 0 ) ≤ C , k ≥ 1 , (2.8)
and the estimate (2.8) implies that

R 2 (u k + v k )(x)(1 + |x| 2 ) dx ≤ C , k ≥ 1 , (2.9)
by classical properties of the Wasserstein distance W 2 , see, e.g., [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF]Eq. (17)]. The compactness of the embedding of

H 1 (R 2 )∩L 1 (R 2 ; |x| 2 dx) in L 2 (R 2
) (which is recalled in Lemma A.1 below) ensures that there are non-negative functions f and g in

H 1 (R 2 ) ∩ L 1 (R 2 , (1 + |x| 2 ) dx
) and a subsequence of (u k , v k ) k≥1 (not relabeled) having the property that

(u k , v k ) → (f, g) in L 2 (R 2 ; R 2 ) and a.e. in R 2 , (u k , v k ) ⇀ (f, g) in H 1 (R 2 ; R 2 ).
(2.10) Furthermore, it readily follows from (2.7), (2.9), and the Dunford-Pettis theorem that (u k ) k≥1 and (v k ) k≥1 are weakly sequentially compact in L 1 (R 2 ). This property, together with (2.10) and the Vitali theorem imply the strong convergence of (u k ) k≥1 and

(v k ) k≥1 in L 1 (R 2 ), namely, (u k , v k ) → (f, g) in L 1 (R 2 ; R 2 ). (2.11)
Since (u k , v k ) k≥1 belongs to K 2 for each k ≥ 1, we conclude from (2.9), (2.10), and (2.11) that (f, g) ∈ K 2 . That (f, g) is a minimizer of F τ in K 2 follows from (2.10), (2.11), the weak lower semicontinuity of E, and the fact that the 2-Wasserstein metric W 2 is lower semicontinuous with respect to the narrow convergence of probability measures in each of its arguments, the latter convergence being guaranteed by (2.11). This completes the proof of the existence of a minimizer to

F τ in K 2 .
Next, when a ≥ b, the uniqueness of the minimizer follows from the convexity of K 2 and W 2 2 and the strict convexity of E.

We finish off the proof by showing that any minimizer (f, g) of F τ in K 2 actually belongs to H 2 (R 2 ; R 2 ), the proof relying on a technique developed in [START_REF] Matthes | A family of nonlinear fourth order equations of gradient flow type[END_REF]. More precisely, denoting the heat semigroup by (G s ) s≥0 which is defined by

(G s h)(x) := 1 4πs R 2 exp - |x -y| 2 4s h(y) dy , (s, x) ∈ [0, ∞) × R 2 , for h ∈ L 1 (R 2 ), classical properties of the heat semigroup ensure that (G s f, G s g) ∈ K 2 for all s ≥ 0 since (f, g) ∈ K 2 . Consequently, F τ (f, g) ≤ F τ (G s f, G s g
) and we deduce that

E(f, g) -E(G s f, G s g) ≤ 1 2τ W 2 2 (G s f, f 0 ) -W 2 2 (f, f 0 ) + B W 2 2 (G s g, g 0 ) -W 2 2 (g, g 0 ) (2.12)
for all s ≥ 0. On the one hand, an explicit computation gives for s > 0

d ds E(G s f, G s g) = R 2 [(A -B) ∇G s f ∂ s (∇G s f ) + B ∇G s (f + g) ∂ s (∇G s (f + g)) +(a -b) G s f ∂ s G s f + b G s (f + g) ∂ s G s (f + g)] dx = -(A -B) ∆G s f 2 2 -B ∆G s (f + g) 2 2 -(a -b) ∇G s f 2 2 -b ∇G s (f + g) 2 2
, which yields, after integration with respect to time,

1 s s 0 (A -B) ∆G σ f 2 2 + B ∆G σ (f + g) 2 2 + (a -b) ∇G σ f 2 2 + b ∇G σ (f + g) 2 2 dσ = E(f, g) -E(G s f, G s g) s .
Since σ → ∆G σ h 2 and σ → ∇G σ h 2 are non-increasing functions for all h ∈ L 1 (R 2 ), we end up with

(A -B) ∆G s f 2 2 + B ∆G s (f + g) 2 2 + (a -b) + ∇G s f 2 2 + b ∇G s (f + g) 2 2 ≤ E(f, g) -E(G s f, G s g) s + (b -a) + ∇f 2 2 (2.13)
for all s > 0. On the other hand, since the heat equation is the gradient flow of the entropy functional H defined by (1.4) for the 2-Wasserstein distance W 2 , see, e.g., [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF][START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF][START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF], it follows from [2, Theorem 11.1.4] that, for all (h, h)

∈ K 2 , 1 2
d ds W 2 2 (G s h, h) + H(G s h) ≤ H( h) for a.e. s ≥ 0. (2.14)
With the choices (h, h) = (f, f 0 ) and (h, h) = (g, g 0 ) in (2.14), we obtain

1 2 d ds W 2 2 (G s f, f 0 ) + B W 2 2 (G s g, g 0 ) ≤ H(f 0 ) -H(G s f ) + B (H(g 0 ) -H(G s g))
for a.e. s ≥ 0. Integrating the above inequality with respect to time and using the time monotonicity of s → H(G s f ) and s → H(G s g) lead us to

1 2s W 2 2 (G s f, f 0 ) -W 2 2 (f, f 0 ) + B W 2 2 (G s g, g 0 ) -W 2 2 (g, g 0 ) ≤ [H(f 0 ) -H(G s f ) + B (H(g 0 ) -H(G s g))] .
(2.15) Combining (2.12), (2.13), and (2.15), we find

(A -B) ∆G s f 2 2 + B ∆G s (f + g) 2 2 + (a -b) + ∇G s f 2 2 + b ∇G s (f + g) 2 2 ≤ 1 τ [H(f 0 ) -H(G s f ) + B (H(g 0 ) -H(G s g))] + (b -a) + ∇f 2 2 (2.16)
for s > 0. Since (f, g) ∈ K 2 , classical properties of the heat semigroup and the functional H entail that (H(G s f ), H(G s g)) converges towards (H(f ), H(g)) as s → 0. This convergence, (2.16), and the positivity of A -B and B readily imply that (∆G s f ) s>0 and (∆G s (f + g)) s>0 are bounded in L 2 (R 2 ). Since they converge towards ∆f and ∆(f + g), respectively, in the sense of distributions as s → 0, we deduce that both ∆f and

∆(f + g) belong to L 2 (R 2 ) and (∆G s f, ∆G s (f + g)) s>0 converges weakly to (∆f, ∆(f + g)) in L 2 (R 2 ; R 2 ) as s → 0. As a consequence, recalling that (f, g) ∈ K 2 , we conclude that f and f + g belong to H 2 (R 2 )
, and so does g. It remains to pass to the limit s → 0 in (2.16) to obtain (2.6).

The Euler-Lagrange equations

We now identify the Euler-Lagrange equation satisfied by the minimizers of the functional

F τ , defined in (2.2), in K 2 . Lemma 3.1. Consider (f 0 , g 0 ) ∈ K 2 and τ > 0. If (f, g) is a minimizer of F τ in K 2 , it satisfies R 2 (f -f 0 )ξ dx + τ R 2 [∆(Af + Bg) div(f ∇ξ) + f ∇(af + bg) • ∇ξ] dx ≤ D 2 ξ ∞ 2 W 2 2 (f, f 0 ) (3.1) and R 2 (g -g 0 )ξ dx + τ R 2 [∆(f + g) div(g∇ξ) + cg ∇(f + g) • ∇ξ] dx ≤ D 2 ξ ∞ 2 W 2 2 (g, g 0 ) (3.2) for ξ ∈ C ∞ 0 (R 2 )
, where D 2 ξ denotes the Hessian matrix of ξ. Proof. By Brenier's theorem [21, Theorem 2.12], there are two measurable functions T : R 2 → R 2 and S : R 2 → R 2 such that (f, g) = (T #f 0 , S#g 0 ) and

W 2 2 (f, f 0 ) = R 2 |x -T (x)| 2 f 0 (x)dx and W 2 2 (g, g 0 ) = R 2 |x -S(x)| 2 g 0 (x)dx . (3.3)
We pick now two test functions η = (η 1 , η 2 ) and ξ = (ξ 1 , ξ 2 ) in C ∞ 0 (R 2 ; R 2 ) and define (T ε , S ε ) := (id +εξ, id +εη) ,

(f ε , g ε ) := (T ε • T )#f 0 , (S ε • S)#g 0 = T ε #f, S ε #g , (3.4)
for each ε ∈ [0, 1], where id is the identity function on R 2 . Clearly, there is ε 0 ∈ (0, 1) small enough (depending on both ξ and η) such that, for ε ∈ [0, ε 0 ], T ε and S ε are C ∞ -diffeomorphisms from R 2 onto R 2 with positive Jacobi determinants

J ξ ε := det(DT ε ) = 1 + ε div(ξ) + ε 2 det(Dξ) , J η ε := det(DS ε ) = 1 + ε div(η) + ε 2 det(Dη) . (3.5)
By (3.4), we have the identities

f ε = f • T -1 ε J ξ ε • T -1 ε and g ε = g • S -1 ε J η ε • S -1 ε , ε ∈ (0, ε 0 ], (3.6)
from which we deduce that f ε and g ε both belong to H 2 (R 2 ) and satisfy

f ε 1 = f 1 = g ε 1 = g 1 = 1. Additionally, we compute f ε 2 2 = R 2 f 2 (x) J ξ ε (x) dx and g ε 2 2 = R 2 g 2 (x) J η ε (x) dx. (3.7)
For further use, we now study the behaviour of (f ε ) ε and (g ε ) ε as ε → 0. To begin with, we notice that, given a test function ϑ ∈ C ∞ 0 (R 2 ), it follows from (3.5) and (3.6) that

lim ε→0 R 2 f ε ϑ dx = lim ε→0 R 2 f (ϑ • T ε ) dx = R 2 f ϑ dx (3.8)
and

lim ε→0 R 2 f ε -f ε ϑ dx = lim ε→0 R 2 f ϑ • T ε -ϑ ε dx = R 2 f ∇ϑ • ξ dx , (3.9)
while (3.5), (3.7), and the Lebesgue dominated convergence theorem entail that

lim ε→0 f ε 2 2 = f 2 2 . (3.10)
Thanks to (3.8) and (3.10) on the one hand and to (3.9) and Lemma A.2 on the other hand, we conclude that

f ε → f in L 2 (R 2 ) and f ε -f ε ⇀ -div(f ξ) in L 2 (R 2 ) . (3.11)
Similarly, we get

g ε → g in L 2 (R 2 ) and g ε -g ε ⇀ -div(gη) in L 2 (R 2 ) . (3.12)
We next turn to the convergence properties of (∇f ε ) ε and (∇g ε ) ε . Differentiating (3.6) and using (3.5), we find

∇f ε • T ε = 1 (J ξ ε ) 2 ∇f + ε V ε (f, ξ) -ε 2 f (J ξ ε ) 3 R ε (ξ) , (3.13)
where

V ε (f, ξ) := (V 1,ε (f, ξ), V 2,ε (f, ξ)), R ε (ξ) := (R 1,ε (ξ), R 2,ε (ξ)), V 1,ε (f, ξ) := ∇f • (∂ 2 ξ 2 , -∂ 1 ξ 2 ) (J ξ ε ) 2 - f ∂ 1 div(ξ) (J ξ ε ) 3 , V 2,ε (f, ξ) := ∇f • (-∂ 2 ξ 1 , ∂ 1 ξ 1 ) (J ξ ε ) 2 - f ∂ 2 div(ξ) (J ξ ε ) 3
, and

R 1,ε (ξ) := ∂ 1 det(Dξ) + ∂ 1 div(ξ) ∂ 2 ξ 2 -∂ 2 div(ξ) ∂ 1 ξ 2 + ε (∂ 1 det(Dξ) ∂ 2 ξ 2 -∂ 2 det(Dξ) ∂ 1 ξ 2 ) , R 2,ε (ξ) := ∂ 2 det(Dξ) + ∂ 2 div(ξ) ∂ 1 ξ 1 -∂ 1 div(ξ) ∂ 2 ξ 1 + ε (∂ 2 det(Dξ) ∂ 1 ξ 1 -∂ 1 det(Dξ) ∂ 2 ξ 1 ) .
Let us now draw several consequences of (3.13): first, we have

∇f ε 2 2 = R 2 |∇f ε • T ε | 2 J ξ ε dx , and, since f ∈ H 1 (R 2 ) and J ξ ε -→ 1 in L ∞ (R 2
) by (3.5), it readily follows from (3.13) and the previous identity that lim

ε→0 ∇f ε 2 2 = ∇f 2 2 . (3.14)
In particular, (∇f

ε ) ε is bounded in L 2 (R 2 ; R 2 ) and, since f ε → f in L 2 (R 2 ) by (3.11), we conclude that (∇f ε ) ε converges weakly towards ∇f in L 2 (R 2 ; R 2 )
. This convergence and (3.14) then guarantee that

∇f ε → ∇f in L 2 (R 2 ; R 2 ) .
(3.15) Next, owing to (3.13), we have

1 ε ∇(f ε -f ) = 1 -J ξ ε •T -1 ε ε ∇f • T -1 ε J ξ ε •T -1 ε 2 + 1 ε ∇f • T -1 ε J ξ ε •T -1 ε -∇f + V ε (f, ξ) • T -1 ε -ε f • T -1 ε J ξ ε •T -1 ε 3 R ε (ξ) • T -1 ε .
The properties of J ξ ε , T ε , f , and the definitions of V ε (f, ξ) and R ε (ξ) readily ensure that the first, third and fourth terms on the right-hand side of the above identity are bounded in L 2 (R 2 ; R 2 ) while the boundedness in L 2 (R 2 ; R 2 ) of the second one follows from Lemma A.

2 since f ∈ H 2 (R 2 ) by Lemma 2.2. Consequently, (∇(f ε -f )/ε) ε is bounded in L 2 (R 2 ; R 2 ) and, recalling that ((f ε -f )/ε) ε converges weakly to -div(f ξ) in L 2 (R 2 ; R 2 ) by (3.11), we conclude that ∇(f ε -f ) ε ⇀ -∇div(f ξ) in L 2 (R 2 ; R 2 ) . (3.16)
Similar computations give

∇g ε → ∇g in L 2 (R 2 ; R 2 ) and ∇(g ε -g) ε ⇀ -∇div(gη) in L 2 (R 2 ; R 2 ) . (3.17)
After this preparation, we can start the proof of (3.1) and (3.2). Recalling that (f ε , g ε ) ∈ K 2 for all ε ∈ (0, ε 0 ], the minimizing property of (f, g) entails that

F τ (f, g) ≤ F τ (f ε , g ε ), that is, 0 ≤ 1 2τ W 2 2 (f ε , f 0 ) -W 2 2 (f, f 0 ) + B W 2 2 (g ε , g 0 ) -W 2 2 (g, g 0 ) + E(f ε , g ε ) -E(f, g). (3.18) Since (T ε • T )#f 0 = f ε by (3.4), we infer from (3.3) that W 2 2 (f ε , f 0 ) ≤ W 2 2 (f, f 0 ) -2ε R (id -T ) • (ξ • T ) f 0 dx + ε 2 R |ξ • T | 2 f 0 dx.
A similar inequality being valid with (g ε , g 0 , g, S ε , S, η) instead of (f ε , f 0 , f, T ε , T, ξ), we conclude that

lim sup ε→0 1 2ε W 2 2 (f ε , f 0 ) -W 2 2 (f, f 0 ) + B W 2 2 (g ε , g 0 ) -W 2 2 (g, g 0 ) ≤ - R 2 (id -T ) • (ξ • T ) f 0 dx + B R 2 (id -S) • (η • S) g 0 dx . (3.19)
We next show that

lim ε→0 1 2ε (a -b) f ε 2 2 + b f ε + g ε 2 2 -(a -b) f 2 2 -b f + g 2 2 = - R 2 a f 2 2 div(ξ) + b g 2 2 div(η) + b f div(gη) + b g div(f ξ) dx . (3.20)
Indeed, we write

(a -b) f ε 2 2 + b f ε + g ε 2 2 -(a -b) f 2 2 -b f + g 2 2 = a I ε 1 + b I ε 2 + b I ε 3 , with I ε 1 := f ε 2 2 -f 2 2 , I ε 2 := g ε 2 2 -g 2 2 , I ε 3 := 2 R 2 (f ε g ε -f g) (x) dx.
It readily follows from (3.11) that

lim ε→0 I ε 1 2ε = lim ε→0 R 2 (f ε + f ) 2 (f ε -f ) ε dx = - R 2 f div(f ξ) dx = - 1 2 R 2 f 2 div(ξ) dx . (3.21)
Similarly, (3.12) guarantees that

lim ε→0 I ε 2 2ε = lim ε→0 R 2 (g ε + g) 2 (g ε -g) ε dx = - 1 2 R 2 g 2 div(η) dx. (3.22)
We next write I ε 3 as

I ε 3 := R 2 [(f ε + f ) (g ε -g) + (g ε + g) (f ε -f )] dx ,
and deduce from (3.11) and (3.12) that Finally, we show that

lim ε→0 I ε 3 2ε = - R 2 [f div(gη) + g div(f ξ)] dx. ( 3
lim ε→0 1 2ε (A -B) ∇f ε 2 2 + B ∇(f ε + g ε ) 2 2 -(A -B) ∇f 2 2 -B ∇(f + g) 2 2 = R 2 [(A∆f + B∆g) div(f ξ) + B∆(f + g) div(gη)] dx. (3.24)
To this end, we write

(A -B) ∇f ε 2 2 + B ∇(f ε + g ε ) 2 2 -(A -B) ∇f 2 2 -B ∇(f + g) 2 2 = A L ε 1 + B L ε 2 + B L ε 3 , with L ε 1 := ∇f ε 2 2 -∇f 2 2 , L ε 2 := ∇g ε 2 2 -∇g 2 2 , L ε 3 := 2 R 2 (∇f ε • ∇g ε -∇f • ∇g) dx.
Thanks to (3.15) and (3.16), we have

lim ε→0 L ε 1 2ε = lim ε→0 R 2 ∇(f + f ε ) 2 • ∇(f ε -f ) ε dx = - R 2 ∇f • ∇div(f ξ) dx = R 2 ∆f div(f ξ) dx , (3.25)
and similarly, by (3.17),

lim ε→0 L ε 2 2ε = lim ε→0 R 2 ∇(g + g ε ) 2 • ∇(g ε -g) ε dx = R 2 ∆g div(gη) dx. (3.26)
Finally,

L ε 3 = R 2 [∇(f ε + f ) • ∇(g ε -g) + ∇(f ε -f ) • ∇(g ε + g)] dx ,
and we infer from (3.15), (3.16), and (3.17) that

lim ε→0 L ε 3 2ε = - R 2 [∇f • ∇div(gη) + ∇div(f ξ) • ∇g] dx = R 2
[∆f div(gη) + ∆g div(f ξ)] dx. We now divide (3.18) by ε and take the limsup as ε → 0, using (3.19), (3.20), and (3.24). The resulting inequality being also valid for (-ξ, -η), we obtain, by choosing successively ξ = 0 and η = 0, that

1 τ R 2 (id -T ) • ξ • T f 0 dx = R 2 [∆(Af + Bg) div(f ξ) + f ∇(af + bg) • ξ] dx , (3.28)
1 τ R 2 (id -S) • η • S g 0 dx = R 2 [∆(f + g) div(gξ) + c f ∇(f + g) • ξ] dx . (3.29) Consider finally Ξ ∈ C ∞ 0 (R 2 ) and take ξ = ∇Ξ in (3.28). Since |Ξ(x) -Ξ(T (x)) -∇Ξ(T (x)) • (x -T (x))| ≤ D 2 Ξ ∞ |x -T (x)| 2 2
for all x ∈ R 2 by the mean value theorem, we find after multiplying the above relation by f 0 and thereafter integrating over R 2 and using (3.3) that

R 2 [Ξ(x) -Ξ(T (x)) -∇Ξ(T (x)) • (x -T (x))] f 0 (x) dx ≤ D 2 Ξ ∞ W 2 2 (f, f 0 ) 2 .
Combining the above inequality with (3.28) gives (3.1). The inequality (3.2) next follows from (3.29) in a similar way.

In order to present the next result, we introduce first some notations. Given a nonnegative and continuous function h and δ > 0, we define the open sets P h δ and P h by P h δ := {x ∈ R 2 : h(x) > δ} and P h := δ>0 P h δ .

Lemma 3.2. Given (f 0 , g 0 ) ∈ K 2 and τ > 0, any minimizer (f, g) of F τ in K 2 is such that Af + Bg ∈ H 3 loc (P f ) and f + g ∈ H 3 loc (P g ). Moreover, the functions j f , w f , j g , and w g defined by

w f := √ f (-∇∆(Af + Bg) + ∇(af + bg)) a.e. in P f , 0 a.e. in R 2 \ P f , , j f := f w f , (3.30)
and

w g := √ g (-∇∆(f + g) + c∇(f + g)) a.e. in P g , 0 a.e. in R 2 \ P g , , j g := √ g w g , (3.31) belong to L 2 (R 2 ; R 2 ) and satisfy R 2 [∆(Af + Bg) div(f ξ) + f ∇(af + bg) • ξ] dx = R 2 j f • ξ dx, (3.32) R 2 [∆(f + g) div(gξ) + c g ∇(f + g) • ξ] dx = R 2 j g • ξ dx, (3.33)
for all ξ ∈ C ∞ 0 (R 2 ; R 2 ). In addition, we have the following estimates: τ w f 2 ≤ W 2 (f, f 0 ) and τ w g 2 ≤ W 2 (g, g 0 ).

(3.34)

Proof. Since H 2 (R 2 ) is embedded in C(R 2
), Lemma 2.2 guarantees that (f, g) ∈ C(R 2 ; R 2 ) so that P f and P g are indeed open subsets of R 2 . Next, recalling (3.28), we use once more the embedding of H 2 (R 2 ) in C(R 2 ) as well as (3.3) to obtain that, for ξ

∈ C ∞ 0 (R 2 ; R 2 ), R 2 [∆(Af + Bg) div(f ξ) + f ∇(af + bg) • ξ] dx ≤ 1 τ R 2 | id -T | 2 f 0 dx 1/2 R 2 |ξ • T | 2 f 0 dx 1/2 ≤ W 2 (f, f 0 ) τ R 2 |ξ| 2 f dx 1/2 ≤ C W 2 (f, f 0 ) τ f 1/2 H 2 ξ 2 .
We may thus extend the functional

ξ -→ R 2 [∆(Af + Bg) div(f ξ) + f ∇(af + bg) • ξ] dx
to a continuous linear functional on L 2 (R 2 ; R 2 ). Consequently, there exists a unique function

j f ∈ L 2 (R 2 ; R 2 ) having the property that R 2 [∆(Af + Bg) div(f ξ) + f ∇(af + bg) • ξ] dx = R 2 j f • ξ dx for all ξ ∈ C ∞ 0 (R 2 ; R 2 ). (3.35)
Since (f, g) ∈ H 2 (R 2 ; R 2 ) by Lemma 2.2, a density argument ensures that the relation (3.35) is actually true for all ξ ∈ H 1 (R 2 ; R 2 ).

Consider now δ > 0 and Ξ ∈ C ∞ 0 (P f δ ; R 2 ). Clearly Ξ/f ∈ H 1 (R 2 ; R 2 ) and we infer from (3.35) with ξ = Ξ/f that

P f δ ∆(Af + Bg) div(Ξ) dx ≤ Ξ L 2 (P f δ ) ∇(af + bg) 2 + Ξ L 2 (P f δ ) j f 2 δ . ( 3
(x) := m 2 χ(mx) for x ∈ R 2 and m ≥ 1. Since H 2 (R 2 ) is embedded in C(R 2 ), we have Y m := 1 m + χ m * f -f 1/2 ∞ -→ 0 as m → ∞ . (3.38) Given ϑ ∈ C ∞ 0 (R 2 ; R 2 ), the vector field ϑ/ √ Y m + χ m * f belongs to C ∞ 0 (R 2 ; R 2 ) too,
= ϑ/ √ Y m + χ m * f , R 2 j f • ϑ √ Y m + χ m * f dx ≤ W 2 (f, f 0 ) τ R 2 |ϑ| 2 f Y m + χ m * f dx 1/2 ≤ W 2 (f, f 0 ) τ f Y m + χ m * f 1/2 ∞ ϑ 2 .
A duality argument then ensures that, for each m ≥ 1,

j f / √ Y m + χ m * f belongs to L 2 (R 2 ; R 2 ) with the estimate j f √ Y m + χ m * f 2 ≤ W 2 (f, f 0 ) τ f Y m + χ m * f 1/2 ∞ . Observing that 0 ≤ f Y m + χ m * f = f -χ m * f + χ m * f Y m + χ m * f ≤ f -χ m * f ∞ Y m + 1 ≤ 1 + Y m ,
we actually have the estimate

j f √ Y m + χ m * f 2 ≤ W 2 (f, f 0 ) τ (1 + Y m ) . (3.39)
Several consequences can be drawn from (3.39): first, since Y m → 0 as m → ∞ by (3.38), the sequence

(j f / √ Y m + χ m * f ) m is bounded in L 2 (R 2 ; R 2
) and there are thus a subsequence of

(j f / √ Y m + χ m * f ) m (not relabeled) and w f ∈ L 2 (R 2 ; R 2 ) such that j f √ Y m + χ m * f ⇀ w f in L 2 (R 2 ; R 2 ) . ( 3

.40)

A simple consequence of (3.38), (3.39), and (3.40) is that

τ w f 2 ≤ W 2 (f, f 0 ) . (3.41)
In addition, since ( √ Y m + χ m * f ) m converges towards √ f uniformly on compact subsets of R 2 , we readily deduce from (3.40) that

j f = f w f a.e. in R 2 .
(3.42)

Next, since f = 0 a.e. in R 2 \ P f , it follows from (3.39) that

R 2 \P f |j f | 2 dx = R 2 \P f |j f | 2 Y m + χ m * f (Y m + χ m * f -f ) dx ≤ j f √ Y m + χ m * f 2 2 (Y m + χ m * f -f ∞ ) ≤ W 2 2 (f, f 0 ) τ 2 (1 + Y m ) 3 Y m -→ m→∞ 0 ,
whence, additionally to (3.37),

j f = 0 a.e. in R 2 \ P f . (3.43)
Finally, owing to (3.40) and (3.43), we have

R 2 \P f |w f | 2 dx = lim m→∞ R 2 \P f w f • j f √ Y m + χ m * f dx = 0 ,
and thus w f = 0 a.e. in R 2 \ P f . This completes the proof of Lemma 3.2 for f . The statements (3.33) and (3.34) for g are proved by similar arguments.

Convergence of the time discretization

We pick now τ > 0 and (f 0 , g 0 ) ∈ K 2 . For each integer n ≥ 1, we define

(f n+1 τ , g n+1 τ ) ∈ K 2 as a solution to the minimization problem inf (u,v)∈K 2 F n τ (u, v) ,
where (f 0 τ , g 0 τ ) := (f 0 , g 0 ) and

F n τ (u, v) := 1 2τ W 2 2 (u, f n τ ) + B W 2 2 (v, g n τ ) + E(u, v) , (u, v) ∈ K 2 . Recall that (f n+1 τ , g n+1 τ
) is well-defined and belongs to H 2 (R 2 ; R 2 ) for all n ≥ 1 by Lemma 2.2. We next let (f τ , g τ ) : [0, ∞) × R 2 → K 2 be the function obtained by the method of piecewise constant interpolation in K 2 as follows: (f τ , g τ )(t) := (f n τ , g n τ ) for all t ∈ [nτ, (n + 1)τ ) and n ∈ N. The next lemma collects estimates which allow us to perform the limit τ → 0 and construct in this way a weak solution of (1.1).

Lemma 4.1.

There is C 3 > 0 depending only on A, B, a, b, c, f 0 , and g 0 such that, for all T ≥ 0 and τ ∈ (0, 1), we have

(i) f τ (T ) 1 = g τ (T ) 1 = 1, (4.1) (ii) ∞ n=0 W 2 2 (f n+1 τ , f n τ ) + W 2 2 (g n+1 τ , g n τ ) ≤ C 3 τ, (4.2) (iii) E(f τ (T ), g τ (T )) ≤ E(f 0 , g 0 ), (4.3)
(iv) R 2 [f τ (T, x) + g τ (T, x)] (1 + |x| 2 ) dx ≤ C 3 (1 + T ), (4.4) (v) max {T,τ } τ ∆f τ (s) 2 2 + ∆g τ (s) 2 2 ds ≤ C 3 (1 + T ), (4.5) (vi) ∞ τ w fτ 2 2 + w gτ 2 2 ds ≤ C 3 , (4.6)
where Proof. The assertion (4.1) follows from the fact that (f n τ , g n τ ) ∈ K 2 for all n ∈ N and τ > 0. We next observe that, since

w fτ := √ f τ (-∇∆(Af τ + Bg τ ) + ∇(af τ + bg τ )) a.
F n τ (f n τ , g n τ ) ≥ F n τ (f n+1 τ , g n+1 τ
) for all n ∈ N, we have

1 2τ W 2 2 (f n+1 τ , f n τ ) + B W 2 2 (f n+1 τ , g n τ ) + E(f n+1 τ , g n+1 τ ) ≤ E(f n τ , g n τ ),
and therefore, for all N ∈ N,

1 2τ N -1 n=0 W 2 2 (f n+1 τ , f n τ ) + B W 2 2 (f n+1 τ , g n τ ) + E(f N τ , g N τ ) ≤ E(f 0 , g 0 ). (4.9)
Recalling that the functional E is bounded from below by Lemma 2.1, we obtain (4.2) after letting N → ∞ in (4.9). Moreover, given T ≥ 0, we choose N ≥ 1 such that T ∈ [N τ, (N + 1)τ ) in (4.9) and arrive at (4.3). Next, the bound (4.4) follows readily from (4.2) and the property (f 0 , g 0 ) ∈ K 2 in a similar manner as (2.9). In order to deduce (4.5), we infer from Lemma 2.2 that, for n ∈ N,

(A -B) ∆f n+1 τ 2 2 + B ∆(f n+1 τ + g n+1 τ ) 2 2 + (a -b) + ∇f n+1 τ 2 2 + b ∇(f n+1 τ + g n+1 τ ) 2 2 ≤ 1 τ H(f n τ ) -H(f n+1 τ ) + B H(g n τ ) -H(g n+1 τ ) + (b -a) + ∇f n+1 τ 2 2 .
Summation from n = 0 to n = N -1 yields

(N +1)τ τ (A -B) ∆f τ (s) 2 2 + B ∆(f τ + g τ )(s) 2 2 +(a -b) + ∇f τ (s) 2 2 + b ∇(f τ + g τ )(s) 2 2 ds ≤ H(f 0 ) -H(f N τ ) + B H(g 0 ) -H(g N τ ) + (b -a) + (N +1)τ τ ∇f τ (s)
(N +1)τ τ (A -B) ∆f τ (s) 2 2 + B ∆(f τ + g τ )(s) 2 2 ds ≤ C H + R 2 f 0 (x) (1 + |x| 2 ) dx + f 0 2 2 + C H + R 2 f N τ (x) (1 + |x| 2 ) dx + B C H + R 2 g 0 (x) (1 + |x| 2 ) dx + g 0 2 2 + C H + R 2 g N τ (x) (1 + |x| 2 ) dx + 4(b -a) + A -B (N +1)τ τ (E(f τ (s), g τ (s)) + C 1 ) ds ≤ C (1 + T ) , for T ∈ [N τ, (N + 1)τ ), hence (4.5). Finally, (4.6) follows from (3.34) and (4.2).
Using uniform estimates from Lemma 4.1, we now establish the time equicontinuity of the family (f τ , g τ ) τ . This step is one of the arguments needed to prove the compactness of (f τ , g τ ) τ . Lemma 4.2. There exists a positive constant C 4 such that, for all t ∈ [0, ∞), s ∈ [0, ∞), and τ ∈ (0, 1) we have

f τ (t) -f τ (s) H -4 + g τ (t) -g τ (s) H -4 ≤ C 4 |t -s| + τ . (4.11)
Proof. Let 0 ≤ s < t with s ∈ [ντ, (ν + 1)τ ), ν ≥ 0, and t ∈ [N τ, (N + 1)τ ), N ≥ ν, be given. By virtue of (3.1), (3.30), and (3.32) we have for

ξ ∈ C ∞ 0 (R 2 ) and n ≥ 1, R 2 (f n τ -f n-1 τ ) ξ dx ≤ τ R 2 j f n τ ∇ξ dx + D 2 ξ ∞ W 2 2 (f n τ , f n-1 τ ) 2 .
Using (3.34), (4.1), and Hölder's inequality, we obtain

R 2 (f n τ -f n-1 τ ) ξ dx ≤ τ f n τ 1/2 1 w f n τ 2 ∇ξ ∞ + D 2 ξ ∞ W 2 2 (f n τ , f n-1 τ ) ≤ ∇ξ ∞ W 2 (f n τ , f n-1 τ ) + D 2 ξ ∞ W 2 2 (f n τ , f n-1 τ ) ,
whence, owing to the continuous embedding of

H 4 (R 2 ) in W 2 ∞ (R 2 ), R 2 (f n τ -f n-1 τ ) ξ dx ≤ C W 2 (f n τ , f n-1 τ ) + W 2 2 (f n τ , f n-1 τ ) ξ H 4 .
Therefore, by (4.2), we have

R 2 (f τ (t) -f τ (s)) ξ dx ≤ N n=ν+1 R 2 (f n τ -f n-1 τ ) ξ dx ≤ C ξ H 4 N n=ν+1 W 2 (f n τ , f n-1 τ ) + W 2 2 (f n τ , f n-1 τ ) ≤ C ξ H 4   √ N -ν N n=ν+1 W 2 2 (f n τ , f n-1 τ ) 1/2 + C 3 τ   ≤ C (N -ν)τ + √ τ ξ H 4 ≤ C √ t -s + τ ξ H 4 ,
which yields (4.11) for f τ . A similar computation based on (3.2), (3.31), and (3.33) gives (4.11) for g τ .

We are now in a position to study the compactness properties of (f τ , g τ ) τ as τ → 0.

Lemma 4.3. There exist nonnegative functions f and

g in C([0, ∞), L 2 (R 2 )
) and a subsequence (τ k ) k≥1 which converges to zero such that, for all t ≥ 0,

(f τ k (t), g τ k (t)) → (f (t), g(t)) in L 2 (R 2 ; R 2 ), (4.12)
(f τ k , g τ k ) → (f, g) in L 2 (0, t; H 1 (R 2 ; R 2 )), (4.13)
and (f (t), g(t)) ∈ K 2 . Moreover, we have (f, g) ∈ L 2 (0, t; H 2 (R 2 ; R 2 )) and

(f τ k , g τ k ) ⇀ (f, g) in L 2 (δ, t; H 2 (R 2 ; R 2 )) (4.14)
for all t > 0 and δ ∈ (0, t).

Proof. On the one hand, we remark that (2.5) together with the estimate (4.3) and Lemma 2.1 imply that (f τ ) τ ∈(0,1) and (g τ ) τ ∈(0,1) are bounded in L ∞ (0, ∞;

H 1 (R 2 ; R 2 )). (4.15)
By interpolation, we have the inequality

h 2 ≤ h 4/5 H 1 h 1/5 H -4 , h ∈ H 1 (R 2 ) ,
which gives, together with (4.11) and (4.15),

f τ (t) -f τ (s) 2 + g τ (t) -g τ (s) 2 ≤ C(|t -s| + τ ) 1/10 (4.16)
for all τ ∈ (0, 1), t ≥ 0, and s ≥ 0. On the other hand, for each t ≥ 0, the sequence (f τ (t), g τ (t)) τ ∈(0,1) lies in a compact subset of L 2 (R 2 ; R 2 ) by (4.4), (4.15), and Lemma A.1. Owing to these two properties, we can invoke [2, Proposition 3.3.1] to conclude that there exists a function (f, g) ∈ C([0, ∞), L 2 (R 2 ; R 2 )) and a subsequence τ k ∈ (0, 1), τ k → 0, such that (4.12) holds true. In addition, we deduce from (4.12), (4.15), and the Lebesgue dominated convergence theorem that

(f τ k , g τ k ) -→ (f, g) in L 2 ((0, T ) × R 2 ) for all T > 0 . (4.17)
We improve now this convergence. Given t ≥ 1 and δ ∈ (0, 1), the estimates (4.4), (4.5), and (4.15) ensure that

t δ f τ (s) 2 H 2 + g τ (s) 2 H 2 ds + sup s∈(δ,t) R 2 (f τ + g τ )(s, x) (1 + |x| 2 ) dx ≤ C (1 + t) . (4.18) By Lemma A.1, H 2 (R 2 ) ∩ L 1 (R 2 , (1 + |x| 2 ) dx) is compactly embedded in H 1 (R 2 ), which in turn is continuously embedded in L 2 (R 2
), and we infer from [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]Lemma 9] that

(f τ k , g τ k ) → (f, g) in L 2 (δ, t; H 1 (R 2 , R 2 )),
which can be improved to (4.13) by using (4.15). Observing next that the right-hand side of (4.18) does not depend on δ, we realize that it follows from (4.18) that, after possibly extracting a subsequence and using a diagonal process, we may assume that (f, g) ∈ L 2 (0, t; H 2 (R 2 ; R 2 )) and that (4.14) holds true. It remains to check that (f (t), g(t)) belongs to K 2 for all t ≥ 0. Owing to (4.12) and (4.15), we readily obtain that f (t) and g(t) are both nonnegative and in H 1 (R 2 ). In addition, (4.12) and (4.18) imply that (f τ k (t), g τ k (t)) k≥1 converges towards (f (t), g(t)) in L 1 (R 2 ; R 2 ) from which we deduce that f (t) 1 = g(t) 1 = 1. Using once more (4.18), this convergence also guarantees that both f (t) and g(t) belong to L 1 (R 2 , (1 + |x| 2 ) dx). Consequently, (f (t), g(t)) ∈ K 2 for all t ≥ 0 and the proof of Lemma 4.3 is complete.

Proof of Theorem 1.1. Let us first check that the functions (f, g) constructed in Lemma 4.3 enjoy the regularity (i) and (ii) stated in Theorem 1.1. The boundedness and integrability properties (i) follow at once from (4.15) and (4.18) by Lemma 4.3. We next use (4.12) to pass to the limit k → ∞ in (4. [START_REF] Matioc | Nonnegative global weak solutions for a degenerate parabolic system modeling thin films driven by capillarity[END_REF]) and obtain

f (t) -f (s) 2 + g(t) -g(s) 2 ≤ C|t -s| 1/10 for all (s, t) ∈ [0, ∞) 2 , (4.19)
which gives the assertion (ii) of Theorem 1.1. We now identify the equations solved by (f, g). For that purpose, we use relations (3.1) and (3.2) to obtain, for N ≥ 1, t ∈ [N τ, (N + 1)τ ), and

ξ ∈ C ∞ 0 (R 2 ), R 2 (f τ (t) -f 0 ) ξ dx + (N +1)τ τ R 2 [∆(Af τ + Bg τ ) div(f τ ∇ξ) + f τ ∇(af τ + bg τ ) • ∇ξ] dx ds ≤ D 2 ξ ∞ 2 N n=1 W 2 2 (f n τ , f n-1 τ ) (4.20) and R 2 (g τ (t) -g 0 ) ξ dx + (N +1)τ τ R 2 [∆(f τ + g τ ) div(g τ ∇ξ) + c g τ ∇(f τ + g τ ) • ∇ξ] dx ds ≤ D 2 ξ ∞ 2 N n=1 W 2 2 (g n τ , g n-1 τ ) . (4.21)
Let now t > 0 be fixed. Before passing to the limit τ → 0 in (4.20) and (4.21), let us point out that, owing to (4.15) and (4.18), we have for all integers ν ≥ 1 and τ > 0 with (ν + 1)τ ≤ t + 1,

(ν+1)τ ντ R 2 [∆(Af τ + Bg τ ) div(f τ ∇ξ) + f τ ∇(af τ + bg τ ) • ∇ξ] dx ds ≤ (ν+1)τ ντ (A ∆f τ 2 + B ∆g τ 2 ) f τ H 1 ξ W 2 ∞ + f τ 2 ∇ξ ∞ (a ∇f τ 2 + b ∇g τ 2 )] ds ≤C(1 + t) √ τ ξ W 2 ∞ (4.22)
and

(ν+1)τ ντ R 2 [∆(f τ + g τ ) div(g τ ∇ξ) + c g τ ∇(f τ + g τ ) • ∇ξ] dx ds ≤ (ν+1)τ ντ ( ∆f τ 2 + ∆g τ 2 ) g τ H 1 ξ W 2 ∞ + c g τ 2 ∇ξ ∞ ( ∇f τ 2 + ∇g τ 2 )] ds ≤C(1 + t) √ τ ξ W 2 ∞ . (4.23)
We fix δ ∈ (0, t). For each k ≥ 1, there are integers N k and ν k such that t ∈ [N k τ k , (N k + 1)τ k ) and δ ∈ [ν k τ k , (ν k + 1)τ k ). In virtue of (4.2), (4.20), and (4.22) we obtain that, for ξ

∈ C ∞ 0 (R 2 ), R 2 (f τ k (t) -f τ k (δ)) ξ dx + t δ R 2 [∆(Af τ k + Bg τ k ) div(f τ k ∇ξ) + f τ k ∇(af τ k + bg τ k ) • ∇ξ] dx ds ≤ D 2 ξ ∞ ∞ p=1 W 2 2 (f p τ k , f p-1 τ k ) + (ν k +1)τ k δ R 2 [∆(Af τ k + Bg τ k ) div(f τ k ∇ξ) + f τ k ∇(af τ k + bg τ k ) • ∇ξ] dx ds + (N k +1)τ k t R 2 [∆(Af τ k + Bg τ k ) div(f τ k ∇ξ) + f τ k ∇(af τ k + bg τ k ) • ∇ξ] dx ds ≤C 3 τ k D 2 ξ ∞ + C(1 + t) √ τ k ξ W 2 ∞ ≤C(1 + t) ξ W 2 ∞ √ τ k . ( 4

.24)

Let us now pass to the limit τ k → 0 in (4.24). We note that the convergences (4.13) and (4.14) guarantee that

∆(Af τ k + Bg τ k ) ∇f τ k ⇀ ∆(Af + Bg) ∇f in L 1 ((δ, t) × R 2 ; R 2 ), ∆(Af τ k + Bg τ k ) f τ k ⇀ ∆(Af + Bg) f in L 1 ((δ, t) × R 2 ), (4.25)
while (4.13) implies

f τ k ∇(af τ k + bg τ k ) -→ f ∇(af + bg) in L 1 ((0, t) × R 2 ; R 2 ). (4.26)
We then let k → ∞ in (4.24) and use (4.12), (4.25), and (4.26) to conclude that

R 2 (f (t) -f (δ)) ξ dx + t δ R 2 [∆(Af + Bg) div(f ∇ξ) + f ∇(af + bg) • ∇ξ] dx ds = 0 for all ξ ∈ C ∞ 0 (R 2 )
. By Lebesgue's dominated convergence theorem and (4.19) we may let δ → 0 and thus obtain the first identity of (1.7). The second identity of (1.7) follows in a similar way, starting from (4.21) and (4.23).

Let us now prove (1.8). We fix t > 0, δ ∈ (0, t) and take k ≥ 1 sufficiently large so that τ k ≤ δ and n k ≥ 1 such that t ∈ [n k τ k , (n k + 1)τ k ). It follows from (4.10) and (4.15) that

t δ D H (f τ k (s), g τ k (s)) ds ≤ (n k +1)τ k τ k D H (f τ k (s), g τ k (s)) + (b -a) ∇f τ k (s) 2 2 ds + (a -b) t δ ∇f τ k (s) 2 2 ds ≤H(f 0 ) -H(f τ k (t)) + B (H(g 0 ) -H(g τ k (t))) + |b -a| δ τ k ∇f τ k (s) 2 2 ds + (n k +1)τ k t ∇f τ k (s) 2 2 ds ≤H(f 0 ) -H(f τ k (t)) + B (H(g 0 ) -H(g τ k (t))) + C |b -a δ . (4.27)
Now, on the one hand, we infer from (4.13) and (4.14) that

lim inf k→∞ t δ D H (f τ k (s), g τ k (s)) ds = lim inf k→∞ t δ D H (f τ k (s), g τ k (s)) + (b -a) ∇f τ k 2 2 ds + lim k→∞ (a -b) t δ ∇f τ k (s) 2 2 ds ≥ t δ D H (f (s), g(s)) + (b -a) ∇f 2 2 ds + (a -b) t δ ∇f (s) 2 2 ds = t δ D H (f (s), g(s)) ds .
On the other hand, it follows from (4.4), (4.12), and (4.15) by classical arguments that

lim k→∞ H(f τ k (t)) + BH(g τ k (t)) = H(f (t)) + BH(g(t)) .
see [START_REF] Ph | A gradient flow approach to a thin film approximation of the Muskat problem[END_REF] for instance. Thanks to these two properties, we can pass to the limit k → ∞ in (4.27) and obtain

t δ D H (f (s), g(s)) ds ≤ H(f 0 ) -H(f (t)) + B (H(g 0 ) -H(g(t))) + C |b -a| δ (4.28)
for all δ < t. By the monotone convergence theorem and the assertion (i) of Theorem 1.1 we may let δ → 0 in (4.28) and end up with (1.8).

In order to obtain the last estimate (1.9), we deduce from (3.34) and (4.9) that, if t > δ > 0 and k is sufficiently large (so that τ k < δ), then

2E(f 0 , g 0 ) ≥ t δ w fτ k 2 2 + B w gτ k 2 2 ds + 2E(f τ k (t), g τ k (t)) , ( 4
) k are bounded in L 2 ((δ, ∞) × R 2 ; R 2
) for all δ > 0. Therefore, after possibly extracting a subsequence and using a diagonal process, we find vector fields V f and

V g in L 2 ((0, ∞) × R 2 ; R 2 ) such that w fτ k , w gτ k ⇀ (V f , V g ) in L 2 ((δ, ∞) × R 2 ; R 2 ) for all δ > 0. (4.30)
Owing to (4.12), (4.15), and (4.30), we can first perform the liminf k → ∞ in (4.29), then take the limit as → 0 with the help of the monotone convergence theorem in the resulting inequality, and thus arrive at

2E(f 0 , g 0 ) ≥ t 0 V f 2 2 + B V g 2 2 ds + 2E(f (t), g(t)
) for all t ≥ 0 . (4.31)

It remains to identify the terms V f and V g . To this end, we remark first that (4.15) ensures that

f τ k ) k and ( √ g τ k ) k are bounded in L ∞ (0, ∞; L 4 (R 2 ; R 2
)), which implies, together with (4.29), that the sequences (j fτ k ) k and (j gτ k ) k defined by

j fτ k = f τ k w fτ k and j gτ k = √ g τ k w gτ k , k ≥ 1, are bounded in L 2 (δ, ∞; L 4/3 (R 2 ; R 2 )) for all δ > 0. Since L 2 (δ, ∞; L 4/3 (R 2
)) is a reflexive space, there are vector fields I f and I g in L 2 (0, ∞; L 4/3 (R 2 ; R 2 )) and a subsequence of (τ k ) k (not relabeled) such that

j fτ k , j gτ k ⇀ (I f , I g ) in L 2 (δ, ∞; L 4/3 (R 2 ; R 2
)) for all δ > 0.

(4.32)

Combining (4.12), (4.30), and (4.32) gives

I f = f V f and I g = √ g V g a.e. in (0, ∞) × R 2 . (4.33) Consider now a test function Ξ ∈ C ∞ 0 ((0, ∞) × R 2 ; R 2 ). For each k ≥ 1 and n ≥ 1, we choose as test function ξ(x) = (n+1)τ k nτ k Ξ(s, x) ds , x ∈ R 2 , in (3.32) for (f n τ k , g n τ k ) and find, since (f τ k , g τ k ) is constant on [nτ k , (n + 1)τ k ): (n+1)τ k nτ k R 2 [∆(Af τ k + Bg τ k ) div(f τ k Ξ) + f τ k ∇(af τ k + bg τ k ) • Ξ] dx ds = (n+1)τ k nτ k R 2 j fτ k • Ξ dx ds .
Summing up the previous identity with respect to n ≥ 1 gives

∞ 0 R 2 [∆(Af τ k + Bg τ k ) div(f τ k Ξ) + f τ k ∇(af τ k + bg τ k ) • Ξ] dx ds = ∞ 0 R 2 j fτ k • Ξ dx ds for k large enough (such that supp(Ξ) ⊂ (τ k , ∞) × R 2 )
. Due to (4.13), (4.14), and (4.32), we can pass to the limit as k → ∞ in the above equality and find

∞ 0 R 2 [∆(Af + Bg) div(f Ξ) + f ∇(af + bg) • Ξ] dx ds = ∞ 0 R 2 I f • Ξ dx ds ,
that is,

I f = -∇(f ∆(Af + Bg)) + ∆(Af + Bg)∇f + f ∇(af + bg) in D ′ ((0, ∞) × R 2 ; R 2 ). ( 4

.34)

A similar argument allows us to deduce from (3.33), (4.13), (4.14), and (4.32) that

I g = -∇(g∆(f + g)) + ∆(f + g)∇g + cg∇(f + g) in D ′ ((0, ∞) × R 2 ; R 2 ). ( 4

Appendix A. Auxiliary results

It is well-known that H 1 (R 2 ) is not compactly embedded in L 2 (R 2 ) due to the non-compactness of R 2 but that compactness can be restored by an additional decay at infinity as in the following lemma:

Lemma A.1. The spaces H 1 (R 2 ) ∩ L 1 (R 2 , |x| 2 dx) and H 2 (R 2 ) ∩ L 1 (R 2 , |x| 2 dx) are compactly embedded in L 2 (R 2 ) and H 1 (R 2 ), respectively. Proof. Let (h k ) k≥1 be a bounded sequence in H 1 (R 2 ) ∩ L 1 (R 2 , |x| 2 dx). Without loss of generality, we may assume that there is a function h ∈ H 1 (R 2 ) such that h k ⇀ h in H 1 (R 2 ). Furthermore, the Rellich theorem guarantees that (h k|D N ) k≥1 is relatively compact in L 2 (D N ) for all integers N ≥ 1, where D N is the open disc centered in zero and of radius N and h k|D N the restriction of h k to D N . We may then extract a subsequence, denoted again by (h k ) k≥1 , such that (h k ) k≥1 (strongly) towards h in L 2 (D N ) for all N ≥ 1.

First, for each N ≥ 1, we have that

h k -h 2 2 ≤ D N |h k (x) -h(x)| 2 dx + R 2 \D N |(h k -h)(x)| 2 dx ≤ D N |(h k -h)(x)| 2 dx + 1 N R 2 \D N |(h k -h)(x)| 2 |x| dx ≤ D N |(h k -h)(x)| 2 dx + 1 N R 2 |(h k -h)(x)| |x| 2 dx 1/2 h k -h 3/2 3 .
Since H 1 (R 2 ) is continuously embedding in L 3 (R 2 ) and the function h belongs to H 1 (R 2 ) and L 1 (R 2 , |x| 2 dx), there exists a constant C such that

h k -h 2 2 ≤ D N |(h k -h)(x)| 2 dx + C N for all k ≥ 1 and N ≥ 1.
Letting first k → ∞ and then N → ∞, we conclude that (h k ) k≥1 converges towards h in L 2 (R 2 ) and thus that H 1 (R 2 ) ∩ L 1 (R and completes the proof.

The next result was used in the identification of the Euler-Lagrange equation for the minimizers of the minimization problem (2.1).

Lemma A.2. Consider h ∈ H 1 (R 2 ), ζ = (ζ 1 , ζ 2 ) ∈ C ∞ 0 (R 2 ; R 2
) and, for ε small enough, define ζ ε := id +εζ, j ε := det(Dζ ε ), and

h ε := (h • ζ -1 ε )/(j ε • ζ -1 ε ). Then, there is ε ζ > 0 such that ((h ε -h)/ε) ε∈(0,ε ζ ) is bounded in L 2 (R 2 ).
Proof. Let us first consider the case h ∈ C ∞ 0 (R 2 ). Then, for ε small enough,

h ε -h 2 2 = R 2 |h ε • ζ ε -h • ζ ε | 2 j ε dx = R 2 h j ε -h • ζ ε 2 j ε dx ≤2 R 2 |h| 2 (1 -j ε ) 2 j ε dx + 2 R 2 |h -h • ζ ε | 2 j ε dx .
Recalling that j ε = 1 + ε divζ + o(ε), we have

1 2 ≤ j ε ≤ 2 and |1 -j ε | ≤ 2 ζ W 1 ∞ ε (A.1)
for ε small enough, and we realize that

2 R 2 |h| 2 (1 -j ε ) 2 j ε dx ≤ 16 ζ 2 W 1 ∞ h 2 2 ε 2 .
Next, using once more (A.1), for sufficiently small ε. Combining the above inequalities gives the claimed boundedness of (h ε -h)/ε in L 2 (R 2 ) for ε small enough. The general case h ∈ H 1 (R 2 ) next follows by a density argument.

2 R 2 |h -h • ζ ε | 2 j ε dx ≤4ε 2
We finally recall some well-known estimates for the functional H defined in (1.4), see, e.g., [

  .23) Combining (3.21), (3.22), and (3.23) gives the claim (3.20).

  25), (3.26), and (3.27) gives the claim (3.24).

  and, by (3.3), (3.28), and (3.35) with the choice ξ

sup k≥1 R 2

 2 |h k (x)| |x| 2 dx ≥ lim k→∞ D N |h k (x)| |x| 2 dx = D N |h(x)| |x| 2 dx, so that h ∈ L 1 (R 2 , |x| 2 dx). Next, we have

ζ 2 ∞ 1 0 R 2 |∇h| 2 j

 2122 ζ sε )(x)| 2 ds dx ≤4ε 2 εs • ζ -1 εs dy ds ≤8ε 2 ζ 2 ∞ ∇h 2 2

  ) be a non-negative function with χ 1 = 1 and set χ m

		.36)
	A duality argument then gives that ∆(Af + Bg) ∈ H 1 (P f δ ) for all δ > 0. Consequently, we get that Af + Bg ∈ H 3 loc (P f ) and together with (3.35) we deduce
	j f = -f ∇∆(Af + Bg) + f ∇(af + bg) a.e. in P f .	(3.37)
	We next prove (3.34), adapting an argument from [18, Proposition 2] and [13, Corollary 2.3]. Let
	χ ∈ C ∞ 0 (R 2

  .29) the functions w fτ k and w gτ k being defined in (4.7) and (4.8), respectively. Since E is bounded from below by Lemma 2.1, we infer from (4.29) that (w fτ k ) k and (w gτ k

  2 , |x| 2 dx) is compactly embedded in L 2 (R 2 ).Consider now a bounded sequence(h k ) k≥1 in H 2 (R 2 ) ∩ L 1 (R 2 , |x| 2 dx).Owing to the previous result, there exist a subsequence, denoted again by (h k ) k≥1 , and a function h∈ H 2 (R 2 ) ∩ L 1 (R 2 , |x| 2 dx) such that (h k ) k≥1 converges towards h strongly in L 2 (R 2 ) and weakly in H 2 (R 2 ). H 2 (R 2 ) ,a simple interpolation argument then gives that (h k ) k≥1 converges towards h strongly in H 1 (R 2 )

	Since	w H 1 ≤ C w	1/2 H 2 w	1/2 2 ,	w ∈

  13, Lemma A.1]. Lemma A.3. Let h be a nonnegative function in L 1 (R 2 , (1+x 2 ) dx)∩L 2 (R 2 ). Then h ln h ∈ L 1 (R 2 )and there exists a positive constant C H such thatR 2 h(x) | ln h(x)| dx ≤ C H +

	R 2	h(x)(1 + |x| 2 ) dx + h 2 2 ,	(A.2)
	H(h) ≥ -C H -

R 2 h(x)(1 + |x| 2 ) dx. (A.3)