On the Muskat flow

Of concern is the motion of two fluids separated by a free interface in a porous medium, where the velocities are given by Darcy's law. We consider the case with and without phase transition. It is shown that the resulting models can be understood as purely geometric evolution laws, where the motion of the separating interface depends in a non-local way on the mean curvature. It turns out that the models are volume preserving and surface area reducing, the latter property giving rise to a Lyapunov function. We show well-posedness of the models, characterize all equilibria, and study the dynamic stability of the equilibria. Lastly, we show that solutions which do not develop singularities exist globally and converge exponentially fast to an equilibrium.


Introduction
The Muskat flow models the evolution of the interface between two fluids in a porous medium and was introduced by Muskat [24] in 1934, see also [25].
Let u i be the velocity, π i the pressure, i the density, and µ i the viscosity of fluid i , respectively. Moreover, let u Γ denote the velocity of Γ = {Γ(t) : t ≥ 0} and V Γ := u Γ · ν Γ the corresponding normal velocity (in the direction of ν Γ ). If there are no sources of mass in the bulk then conservation of mass is given by the continuity equation If there is no surface mass on Γ(t), we also have the jump condition We note that j Γ is well-defined, as (1.1) shows. If j Γ ≡ 0 then and in this case, the interface Γ(t) is advected with the velocity field u. On the other hand, if j Γ ≡ 0, phase transition occurs, and the normal velocity can then be expressed as ]. In this case we will always assume that 1 = 2 . In the following, we will only consider the completely incompressible case where i > 0 is constant.
Modeling flows in porous media often relies on Darcy's law, which reads where κ > 0 is the permeability of the porous medium; to shorten notation, we set k i = κ/µ i , i = 1, 2. If no phase transition takes place we obtain from Darcy's law (1.4) and the normal velocity is then given by (1.5) By (1.4), the right hand side of (1.5) does not depend on the phases, and hence the expression for V Γ is unambiguous. In case of phase transition, the normal velocity is given by Finally, we assume that the capillary pressure π c = [[π]] is given by where H Γ = −div Γ ν Γ denotes the (n − 1)-fold mean curvature, that is, the sum of the principal curvatures of Γ(t), and σ > 0 is the surface tension. Here the convention is that H Γ = −(n − 1)/R for a sphere of radius R in R n .
The resulting problem in the case without phase transition is the well-known Muskat problem, or Muskat flow, which is given by ∆π = 0 in Ω \ Γ(t), (1.8) If j Γ = 0 and 1 = 2 , we obtain the Muskat flow with phase transition ∆π = 0 in Ω \ Γ(t), (1.9) see [27, Chapter 1] for a derivation. In fact, in [27] the more general situation where the motion of the fluids is governed by the Navier-Stokes equations is also considered.
For later use we note that the scaled function p = π/ , with π a solution of (1.9), satisfies the equivalent problem ∆p = 0 in Ω \ Γ(t), H Γ on Γ(t), (1.10) In more generality, one can also consider the case where κ i depends on the pressure π i . A variant of Darcy's law is Forchheimer's law, which reads where the function g is strictly positive and s → sg(s) is strictly increasing. Solving this equation for u i one obtains where k i is strictly positive and satisfies k i (p, s) + 2s∂ 2 k i (p, s) > 0 for p ∈ R and s ≥ 0. These conditions ensure strong ellipticity of the second-oder differential operator −div(k i (π i , |∇π i | 2 )∇ ·) in Ω i . In case of non-constant densities = (π), the first line in (1.8) and (1.9) ought to be replaced by ∂ t (π) − div (π)k(π, |∇π| 2 )∇π = 0, π(0) = π 0 .
The resulting model is known as the Verigin problem (with phase transition in case j Γ = 0.) This problem is studied in [28].
It will be shown in Section 3 that problems (1.8) and (1.9) can be cast as a geometric evolution equation where one aims to find a (sufficiently smooth) family of hypersurfaces Γ(t) ⊂ Ω which enclose a domain Ω 1 (t). Here G Γ : W Suppose that the disperse region Ω 1 consists of m ≥ 1 connected components Ω 1,j , that is, Ω 1 = ∪ m j=1 Ω 1,j , while Ω 2 is connected, see Figure 1. Let E denote the set of equilibria for (1.8) and (1.9).  It is interesting to note that the Mullins-Sekerka problem, given by enjoys the same geometric properties as the Muskat flow with phase transition (1.9), see [20,27]. Problems (1.8), (1.9), and (1.12) are all of order 3, and their principal linearizations have equivalent symbols. Finally, we note that the Mullins-Sekerka problem (1.12) has the same set of equilibria as problem (1.9), with analogous stability properties. Problem (1.12) is also known as the quasi-stationary Stefan problem with surface tension and it describes the motion of a material with phase transition, with θ the temperature and d i the respective (constant) diffusion coefficients. We refer again to the monograph [27] for a comprehensive discussion of the physical background.
The Muskat problem has recently received considerable attention. In the case of σ > 0, the first result on the existence of classical solutions in two dimensions was obtained by Hong, Tao and Yi [22]. Regarding the stability of equilibria, Friedmann and Tao [21] proved stability of a circular steady-state in case that Ω 2 is unbounded. The authors of [13] state that the equilibrium is in general not asymptotically stable.
Escher and Matioc [18] considered the Muskat problem in a horizontally periodic geometry with surface tension and gravity included. Existence and uniqueness of classical solutions is obtained and the authors establish exponential stability of certain flat equilibria. Using bifurcation theory they also identify finger shaped steady-states which are all unstable. These results were later refined and extended in Ehrnström, Escher, Matioc, Walker [16,17,19]. Bazaliy and Vasylyeva [2] first observed a waiting time behavior for the two-dimensional Muskat problem with a non-regular initial surface in the presence of surface tension.
The Muskat problem with phase transition (1.9) has been introduced for the first time in [27].
Throughout this paper, we use the notation B X (x, r) for a ball or radius r and center x, with X a normed vector space. For two given normed vector spaces X and Y , B(X, Y ) denotes the space of all bounded linear operators from X into Y , equipped with the uniform operator norm.

Elliptic transmission problems
In this section we consider an elliptic transmission problem which turns out to be important for the analysis of the Muskat flow (1.8).
Suppose that Ω ⊂ R n is a bounded domain with C 2 -boundary, consisting of two parts Ω 1 and Ω 2 , as depicted in Figure 1. Moreover, suppose that Γ = ∂Ω 1 is The following elliptic transmission problem, whose formulation is more general than actually needed for this paper, is also of independent interest.
Proposition 2.1. Let 1 < p < ∞. Then there exists ω 0 ∈ R such that the transmission problem (2.1) has for each ω > ω 0 and each Proof. Here we give a sketch of the proof, and refer to [27,Chapter 6] for more details.
(a) We first consider the case with constant coefficients a 1 , a 2 , flat interface Γ = with ν = e n the outer unit normal of Ω 1 .
To obtain solvability of the problem in the right regularity class, we transform the problem to the half-space case as follows. Set for (x, y) ∈ R n−1 × (0, ∞), and consider the problem where the subscripts 1, 2 refer to the coefficients in the lower resp. upper half-space. Problem (2.3) is strongly elliptic and satisfies the Lopatinskii-Shapiro condition for the half space. By well-known results for elliptic systems, see for instance [27, Section 6.3], this problem is uniquely solvable in the right class, hence the transmission problem (2.2) has this property as well. This proves Proposition 2.1 for the constant coefficient case with flat interface.
(b) By perturbation, the result for the flat interface with constant coefficients remains valid for variable coefficients with small deviation from constant ones. By another perturbation argument, a proper coordinate transformation transfers the result to the case of a bent interface. The localization technique finally yields the result for the case of general domains and general coefficients, see for instance [27] Section 6.3 for more details.
The transmission problem ∆u = 0 in Ω \ Γ, is closely related to the Muskat problem (1.8). As in Section 1, k i = k| Ωi is assumed to be constant for i = 1, 2. We have the following result on solvability. Proof. By Proposition 2.1 we know that problem (2.4), with the first line replaced by provided ω 1 is sufficiently large. In addition, one readily verifies that u 1 ∈ L p,0 (Ω). Next we show that the problem −k∆ũ = ω 1 u 1 in Ω \ Γ, has a unique solutionũ ∈ W 2 p (Ω \ Γ) ∩ L p,0 (Ω). In order to see this, let X = L p,0 (Ω) and consider the linear operator L : D(L) ⊂ X → X given by Then L has compact resolvent and therefore, its spectrum consists only of eigenvalues of finite algebraic multiplicity which, in addition, do not depend on p. By a standard energy argument we obtain σ(L) ⊂ R + . The fact that we restrict ourselves to functions with mean zero implies that 0 lies in the resolvent set of L. Therefore, (2.5) has a unique solutionũ ∈ W 2 p (Ω\Γ)∩L p,0 (Ω). It is now clear that the function u = u 1 +ũ satisfies the assertions of the proposition.

Volume, area, and equilibria
In this section we show that the Muskat problems (1.8) and (1.9) enjoy some important geometric properties, namely conservation of volume and decrease of surface area. Moreover, we characterize all the equilibria. We start by showing that both problems (1.8) and (1.10) can be rewritten as a geometric evolution equation for the motion of Γ(t). Here is linear and satisfies This can be seen as follows. Given h ∈ W It is now clear that (3.1) is equivalent to (1.8). For g, h ∈ W 3/2 2 (Γ) let p(g), p(h) ∈ W 2 2 (Ω\Γ)∩L 2,0 (Ω) be the corresponding solutions of (3.3). Then one readily verifies that showing that G Γ satisfies (3.2). For the Muskat problem with surface tension (1.10) we proceed as follows. Given h ∈ W (3.6) we see that the Muskat problem (1.10) can be rewritten as (3.1). For g, h ∈ W 3/2 2 (Γ), let p i (g), p i (h) ∈ W 2 2 (Ω i ) be the corresponding solutions of (3.5) and (3.6), respectively. Then one verifies that and this shows that (3.2) also holds for (1.10).
Let Ω 1,j , j = 1, . . . , m, denote the components of Ω 1 and Γ j their boundaries, and let Γ := We are now ready for the proof of Theorem 1.1(a)-(b): (a) Let |Ω 1 (t)| denote the volume of Ω 1 (t). By the change of volume formula, see for instance [27, Section 2.5], and (3.9) we obtain For problem (1.8) we obtain by (3.8) showing that |Γ| is decreasing, and hence is a Lyapunov function, for (1. showing that p is constant on Ω 2 and on the connected components Ω 1,j of Ω 1 . Therefore, h is constant on Γ j , that is, h = m j=1 a j e j with some real numbers a j . This implies that H Γ is constant on each component Γ j of Γ, and by Alexandrov's characterization of compact closed hypersurfaces with constant mean curvature, Γ is the union of disjoint spheres, which may all have different radii. This, in turn, also yields that the equilibria for (1.8) consist of disjoint spheres of arbitrary radii. One shows that E, the set of all equilibria, is a smooth manifold of dimenion m(n + 1), see for instance [20,27].
For the Muskat problem with phase transition we proceed analogously: let p i , i = 1, 2, be the solution of (3.5) and (3.6), respectively. Then (3.7) implies that p i is constant on the connected components of Ω. The condition [[p]] = 0 in turn shows that p ≡ c on Ω, and this implies that the mean curvature H Γ is constant all over Γ. Consequently, Γ is the disjoint union of spheres of the same radius and E has dimension mn + 1.

Well-posedness
In this section, we show that the evolution equation (3.1) admits a unique solution which instantaneously regularizes, provided Γ 0 ∈ W s p with s > 2 + (n − 1)/p. In order to establish this result, we use the common approach of transforming problems (1.8) and (1.10), or equivalently problem (3.1), to a domain with a fixed interface Σ, where Γ(t) is parameterized over Σ by means of a height function h(t). For this we rely on the Hanzawa transform, see for instance [27,Section 1.3.2].
We assume, as before, that Ω ⊂ R n is a bounded domain with boundary ∂Ω of class C 2 , and that Γ ⊂ Ω is a hypersurface of class C 2 , i.e., a C 2 -manifold which is the boundary of a bounded domain Ω 1 ⊂ Ω. As above, we set Ω 2 = Ω\Ω 1 , see again Figure 1. If follows from the results in [27,Section 2.3.4], see also [26], that Γ can be approximated by a real analytic hypersurafce Σ, in the sense that the Hausdorff distance of the second order normal bundles is as small as we please. More precisely, given η > 0, there exists an analytic hypersurface Σ such that d H (N 2 Σ, N 2 Γ) ≤ η. If η > 0 is small enough, then Σ bounds a domain Ω Σ 1 with Ω Σ 1 ⊂ Ω and then we set Ω Σ 2 = Ω \ Ω Σ 1 ⊂ Ω. In the sequel we will freely use the results from [27,Chapter 2]. In particular, we know that the hypersurface Σ admits a tubular neighborhood, which means that there is a 0 > 0 such that the map Λ : Σ × (−a 0 , a 0 ) → R n Λ(p, r) := p + rν Σ (p) is a diffeomorphism from Σ × (−a 0 , a 0 ) onto im(Λ), the image of Λ. The inverse of this map is conveniently decomposed as Here Π Σ (x) means the metric projection of x onto Σ and d Σ (x) the signed distance from x to Σ; so |d Σ (x)| = dist(x, Σ) and d Σ (x) < 0 if and only if x ∈ Ω Σ 1 . In particular we have im(Λ) = {x ∈ R n : dist(x, Σ) < a 0 }. The maximal number a 0 is given by the radius r Σ > 0, defined as the largest number r such the exterior and interior ball conditions for Σ in Ω holds.
If dist(Γ, Σ) is small enough, we may use the map Λ to parameterize the unknown free boundary Γ(t) over Σ by means of a height function h(t) via for small t ≥ 0, at least. We then extend this diffeomorphism to all ofΩ by means of a Hanzawa transform. With the Weingarten tensor L Σ and the surface gradient ∇ Σ we further have The transformed problem then reads with µ ∈ (1/3 + (n + 3)/3p, 1]. Here we note that this choice of µ implies the embedding X γ,µ → C 2 (Σ), showing that the mean curvature H Γ (h) is well-defined.
Theorem 4.1. Let p ∈ (1, ∞), and let the spaces X 0 , X 1 , and X γ,µ be defined as in (4.2). Then (3.1) is locally well-posed in the sense that the transformed problem (4.1) is locally well-posed for initial values h 0 ∈ X γ,µ which are small in the topology of C 1 (Σ). Furthermore, the map t → Γ(t) is real analytic.
Proof. We want to rewrite (4.1) as a quasilinear evolution equation where h 0 is small in C 1 (Σ). We recall the representation of the curvature H Γ from [27, Section 2.2.5], which reads , where c 0 and c 1 are real analytic functions, c 0 (0, 0) = I, c 1 (0, 0) = H Σ , and −H Γ is strongly elliptic if h is small in C 1 (Σ). Next one shows that the map is real analytic, provided h is small with respect to the topology of C 1 (Σ). Furthermore, we write β(h) −1 G Γ (h) = G Σ (h), resulting in the problem Here we note that G Σ is a linear pseudo-differential operator of order 1 on Σ for both Muskat problems (1.8) and (1.10). We use the decomposition By the techniques developed in [27,Section 9.5], it is not difficult to show that is real analytic, provided r > 0 is small enough. Key for this is the embedding X γ,µ → C 2 (Σ) which is ensured by the choice of µ. It remains to show that A(h) has the property of L p -maximal regularity.
In order to see this, we note that where ∆ Σ is the Laplace-Beltrami operator on Σ. It follows from Corollaries 6.6.5 and 6.7.4 in [27] that the operator −A(h) with domain D(A(h)) = X 1 has L pmaximal regularity in X 0 for both problems (1.8) and (1.9) for each h ∈ B Xγ,µ (0, r), provided r is sufficiently small. Therefore, Theorems 5.1.1 and 5.2.1 in [27] apply to obtain local well-posedness as well as analyticity in time. For analyticity in space we may follow the arguments presented in [27,Section 9.4].

Stability of equilibria
Recall that the equilibria of (1.8) and (1.9) consist of finitely many spheres Σ j := S(x j , R j ), 1 ≤ j ≤ m. Given such an equilibrium Γ * = m j=1 Σ j , we choose Σ = Γ * as the reference hypersurface. The linearization of the transformed problem then reads with R j the radius of the sphere Σ j , and ∆ Σj the Laplace-Beltrami operator of Σ j . This follows from the fact that the Fréchet derivative of G Σ (h)H Γ (h) at h = 0 (in the direction of g) can be evaluated by as H Σ = H Γ (0) is constant on equilibria, and G Σ (εg)e = 0. As the operator −G Σ A Σ has maximal regularity, we may apply the stability results from [27, Chapter 5], once we have shown that 0 is normally stable or normally hyperbolic for (4.3).
Before showing the latter we recall the pertinent definitions. Let L := σG Σ A Σ be the linearization of −σG Σ (h)H Γ (h) at the equilibrium h = 0.
We are ready to prove the following important result.
Proof. It follows from our previous considerations that the set of of equilibria form a smooth manifold. Next we note that G Σ A Σ has compact resolvent by boundedness of Ω, so we only need to consider its eigenvalues.
(a) We begin with eigenvalue 0. So let G Σ A Σ h = 0. Then A Σ h belongs to the kernel of G Σ , which implies by (3.9) that A Σ h = ae in case (1.9), and A Σ h = m j=1 a j e j in case (1.8), see (3.8). Therefore, h = h 0 − (R 2 /(n − 1))ae for (1.9), and h = h 0 − m j=1 (R 2 j /(n − 1))a j e j in case of (1.8), where h 0 ∈ N(A Σ ). As dim N(A Σ ) = mn, we conclude that the dimension of the kernel N(G Σ A Σ ) equals the dimension of the manifold E.
(b) To see that the eigenvalue 0 is semi-simple for G Σ A Σ , suppose (G Σ A Σ ) 2 h = 0. Then for (1.8) a j e j , for some h 0 ∈ N(A Σ ), a j ∈ C.
Multiplying this relation with e l in L 2 (Σ) we obtain a j = 0 for all j, as G Σ is selfadjoint and G Σ e j = 0. As A Σ is also selfadjoint, multiplying with A Σ h, we obtain (G Σ A Σ h|A Σ h) L2(Σ) = 0, hence G Σ A Σ h = 0. The argument for (1.9) is similar. Consequently, 0 is semi-simple for G Σ A Σ .
As G Σ and A Σ are selfadjoint, this identity implies that λ must be real, hence the spectrum of G Σ A Σ is real.
We consider now the case (1.8); then (h|e j ) Σ = 0 for all j = 1, . . . , m. Suppose λ > 0. As G Σ is positive semi-definite and A Σ is so on the orthogonal complement of span{e j } m j=1 we see that (h|A Σ h) = 0. This implies A Σ h = 0 and then h = 0 as λ > 0. Therefore, there are no nonzero eigenvalues with nonnegative real part, hence in this case 0 is normally stable.
In case (1.9), we only obtain (h|e) Σ = 0. As A Σ is positive semi-definite on functions with mean zero if and only if Σ is connected, we may conclude normal stability, provided Σ is connected.
• Regularity: the norm of Γ(t) in W s p may become unbounded as t → t + (Γ 0 ); • Geometry: the topology of the interface Γ(t) may change, or the interface may touch the boundary of Ω, or part of it may shrink to points. We say that the solution Γ(t) satisfies auniform ball condition, if there is a number r > 0 such that for each t ∈ J 0 := [0, t + (Γ 0 )) and each p ∈ Γ(t) there are balls B(x i , r) ⊂ Ω i , i = 1, 2, such thatB(x i , r) ∩ Γ(t) = {p}. The main result of this section reads as follows.
Theorem 6.1. Let Γ(t) be a solution of the geometric evolution equation (3.1) on its maximal time interval J(Γ 0 ). Assume furthermore that (i) |Γ(t)| W s p ≤ M < ∞ for all t ∈ J(Γ 0 ), and (ii) Γ(t) satisfies a uniform ball condition. Then J(Γ 0 ) = R + , i.e., the solution exists globally, and Γ(t) converges in SM s to an equilibrium Γ ∞ ∈ E at an exponential rate. The converse is also true: if a global solution converges in SM s to an equilibrium, then (i) and (ii) are valid.
Proof. The proof relies on [23,Theorem 4.3] and follows the same lines as that of Theorem 5.2 in [23]; see also [27], Theorems 5.7.2 and 11.4.1.