## Abstract

In this note, we show that there exist solutions of the Muskat problem that shift stability regimes: they start unstable, then become stable and finally return to the unstable regime. We also exhibit numerical evidence of solutions with medium-sized norm of the derivative of the initial condition that develop a turning singularity.

## 1. Introduction

In this paper, we study two incompressible fluids with the same viscosity but different densities, *ρ*^{+} and *ρ*^{−}, evolving in a two-dimensional porous medium with constant permeability *κ*. The velocity *v* is determined by Darcy's law
1.1
where *p* is the pressure, *μ*>0 viscosity, and *g*>0 gravitational acceleration. In addition, *v* is incompressible
1.2
By rescaling properly, we can assume *μ*=*g*=1. The fluids also satisfy the conservation of mass equation
1.3

This is known as the Muskat problem [1]. We denote by *Ω*^{+} the region occupied by the fluid with density *ρ*^{+} and by *Ω*^{−} the region occupied by the fluid with density *ρ*^{−}≠*ρ*^{+}. The point belongs to *Ω*^{+}, whereas the point belongs to *Ω*^{−}. All quantities with superindex ± will refer to *Ω*^{±}, respectively. The interface between both fluids at any time *t* is a curve *z*(*α*,*t*). We will work in the setting of flat at infinity interfaces, although the results can be extended to the horizontally periodic case.

A quantity that will play a major role in this paper is the Rayleigh–Taylor (RT) condition, which is defined as
where we use the convention (*u*,*v*)^{⊥}=(−*v*,*u*). If RT(*α*,*t*)>0 for all , we will say that the curve is in the Rayleigh–Taylor stable regime at time *t*, and if *RT*(*α*,*t*)≤0 for some , we will say that the curve is in the Rayleigh–Taylor unstable regime.

One can rewrite the system (1.1)–(1.3) in terms of the curve *z*(*α*,*t*), obtaining
1.4
A simple calculation of the Rayleigh–Taylor condition in terms of *z*(*α*,*t*) yields
When the interface is a graph, parametrized as *z*(*α*,*t*)=(*α*,*f*(*α*,*t*)), equation (1.4) becomes
1.5
and the Rayleigh–Taylor condition simplifies to
The curve is now in the Rayleigh–Taylor stable regime whenever *ρ*^{+}<*ρ*^{−}, i.e. the denser fluid is at the bottom.

The Muskat problem has been studied in many works. A proof of local existence of classical solutions in the Rayleigh–Taylor stable regime and ill-posedness in the unstable regime appears in [2]. A maximum principle for can be found in [3]. Moreover, the authors showed in [3] that if , then for all *t*>0. Further work has shown existence of finite time turning (i.e. the curve ceases to be a graph in finite time and the Rayleigh–Taylor condition changes sign) [4]. The precise result is the following.

### Theorem 1.1

*There exists a non-empty open set of initial data in H*^{4} *with Rayleigh–Taylor strictly positive (i.e. RT(α,0)>0 for all* *) such that the solutions of (1.4) have RT(α,t)<0 for some finite time t>0 and all α in some non-empty open interval.*

After this shift of regime, the curve may lose regularity. This was proved in [5]. More general models which take into account finite depth or non-constant permeability that also exhibit turning were studied in [6,7], where the estimates in the latter one were carried out by rigorous computer-assisted techniques, as opposed to the traditional pencil and paper ones from the former. All these results are local in time, and the techniques employed to prove them rely on a local analysis of the equations; therefore, it is not possible to conclude global properties of the solutions from them. Concerning global existence, the first proof for small initial data was carried out in [8] in the case where the fluids have different viscosities and the same densities (see also [2] for the setting of this paper: different densities and the same viscosities). Global existence for medium-sized initial data was established in [9,10]. In the case where surface tension is taken into account, global existence was shown in [11,12]. Global existence for the confined case was treated in [13]. Recently, in [14], a new framework was used to study global existence. The following theorem was proved in [9]:

### Theorem 1.2

*Suppose that* *and* *. Then there exists a global in time weak solution of (1.5) that satisfies*
*In particular, f is Lipschitz continuous.*

The question of whether the constant 1 (which is independent of the physical parameters of the system) in the bound of the derivative of the initial data is sharp or not is still open. The proof of theorem 1.1 constructs initial conditions that have parameters that need to be taken big enough to ensure turning. Taking both these two facts into account, our initial motivation was to quantify and bridge the gap, finding or approximating the threshold *C* for the norm of ∂_{x}*f*_{0} such that initial data for which exist globally, and for each *ϵ*>0 there exist solutions with which develop turning singularities.

This note is organized as follows: in §2, we describe the numerical experiments that became the motivation for our main theorem, which we prove in §3. Finally, in §4 we discuss open problems and future work.

## 2. Numerical results

In this section, we describe numerical experiments which led to our main result, theorem 3.3.

Our motivation was to find initial data, as small as possible, which develop turning in finite time. For simplicity, we performed simulations in the horizontally periodic scenario. The evolution equation for the interface reads in that case
2.1
Picking different initial conditions and evolving them until they either develop turning or flatten out seems like looking for a needle in a haystack. Instead, we took data which we were sure that would turn (the interface has a vertical tangent at a point and the velocity is pointing in the right direction) and ran the equation *backwards* in time for a short time to find our desired initial condition.

At the linear level and backwards in time, equation (2.1) behaves like a backwards heat equation, which is ill posed if the lighter fluid is on top of the denser one. To overcome this difficulty, we use the following heuristic: we perform very small steps backwards in time and at every step, we smooth the function by eliminating all frequencies whose components are below a given threshold. The reason behind this heuristic is that the family of solutions which turnover is an open set, and therefore small perturbations (the regularized versions) of the solution should remain in it. We remark that the backwards evolution is not done with the purpose of finding a numerical solution of the equation, but only to gain intuition about what initial condition to choose. Once we find a suitable candidate, we check its validity by running the equation forward with the candidate.

The smoothing threshold was set to 10^{−6}. The time integration was done using a Runge–Kutta 45 scheme. The derivatives were calculated using a spectral method with a cutoff filter given by
In order to evaluate the singular integrals, we used an alternating quadrature rule [15,16]
where *h*=2*π*/*N* and *N* is the total number of nodes. We chose the condition at time *t*=0 to be
and *N*=2048 nodes (figure 1).

The smoothing was done after every *Δt*=4×10^{−5} and the simulation ran until *t*_{f}=−4.92×10^{−2}. The obtained evolution of *z*(*α*,*t*) is depicted in figure 2.

Finally, we computed the evolution of equation (2.1) forward in time taking as initial condition *z*(*α*,−4.92×10^{−2}), which has two vertical tangents at approximately (±3.795×10^{−3},±1.268). The initial data are now in the stable regime and the equation is well posed. We can use the same integration scheme (both in time and space) with no smoothing. We can see that the curve comes from the unstable regime, then moves to the stable one, and finally returns to the unstable regime. The whole simulation is carried out in the stable regime.

## 3. Statement and proof of the main result

Our main result, theorem 3.3, is motivated by the above numerics (although the mechanics of the constructed solutions will be somewhat different). This section is devoted to its proof, which will rest on the following two lemmas.

### Lemma 3.1

*There exists an analytic, odd, asymptotically flat initial condition* *z*(*α*,0) *whose analytic extension is* *H*^{4} *on the boundary of its strip of analyticity, and also* *α*_{1}>0 *such that the following hold. We have* ∂_{α}*z*_{1}(*α*,0)=*z*_{2}(*α*,0)=0< ∂_{α}*z*_{2}(*α*,0) *for* *α*∈{0,±*α*_{1}}, *while* ∂_{α}*z*_{1}(*α*,0)>0 *for all other* (*in particular*, *z*(*α*,0) *is a graph of a function of* *x* *with three vertical tangents*), *such that the corresponding solution* *z*(*α*,*t*) *of* (1.4) *satisfies* sgn(*t*)∂_{α}*z*_{1}(±*α*_{1},*t*)>0 *and* sgn(*t*)∂_{α}*z*_{1}(0,*t*)<0 *for all times* *t* *sufficiently close to* 0.

By asymptotically flat, we mean *z*(*α*,0)≈(*α*,0) and ∂_{α}*z*(*α*,0)≈(1,0) for |*α*|≫1. Therefore, *z*(*α*,*t*) will be a graph of a function of *x* everywhere except near *α*=0 for small *t*>0, and everywhere except near *α*=±*α*_{1} for small *t*<0.

### Proof.

Since in the following we will mostly consider *t*=0, let us denote *z*(*α*)= *z*(*α*,0). Following the arguments of [4], we find that if *α*_{0} is any point with *z*_{1}′(*α*_{0})=*z*_{1}′′(*α*_{0})=*z*_{2}(*α*_{0})=0, then the corresponding velocity *v* at time *t*=0 satisfies
We obviously have the following:

— if ∂

_{α}*v*_{1}(*z*(*α*_{0}))>0, then sgn(*t*)∂_{α}*z*(*α*,*t*)>0 for (*α*,*t*) close to (*α*_{0},0). In particular, the curve will be in the stable regime near*α*_{0}for small*t*>0; and— if ∂

_{α}*v*_{1}(*z*(*α*_{0}))<0, then sgn(*t*)∂_{α}*z*(*α*,*t*)<0 for (*α*,*t*) close to (*α*_{0},0). In particular, the curve will be in the unstable regime (i.e. it will turnover) near*α*_{0}for small*t*>0.

Our *z*(*α*) will consist of three building blocks: centre and two tails. First let
and
(note that both are odd). We will now consider curves *z*^{R}(*α*) of the following form:
where *R*>9 will be fixed later. A sketch of the initial condition *z*^{R}(*α*), which satisfies all hypotheses of the lemma (with *α*_{1}=*R*) except regularity, appears in figure 3. We will now show that for small *t*>0, the corresponding solution of (1.4) will turnover near *α*=0 but not near *α*=±*R*. To do so, we will split the integrals in the formula for ∂_{α}*v*_{1}(*z*^{R}(*α*_{0})) (*α*_{0}=0,±*R*) into centre–centre, tail–tail, centre–tail and tail–centre contributions (note also that ). We will compute the first two and only estimate the last two since they will be made arbitrarily small by choosing *R* large enough.

We start with centre–centre. Since *z*^{c} is odd, the quantity we are interested in simplifies to
where
The last three integrals can be explicitly calculated as follows:
The first one can also be calculated explicitly, with slightly more effort. For the sake of simplicity, we do not present here the full symbolic integral, only an approximation of the final result
Adding all the contributions, we obtain

Next, the tail–centre term (i.e. the contribution of the tails to ∂_{α}*v*_{1}(*z*^{R}(0))) is
The explicit expression for *z*^{t} yields the easy bound
Hence, by choosing *R* large enough we can ensure that
3.1

Let us now consider the tail–tail contribution. Because of symmetry it suffices to consider *α*_{0}=*R*. Then we have
where (after using the explicit expression for *z*^{t})
The first and last integral can again easily be evaluated explicitly
and
whereas since its integrand is positive. We can thus conclude that .

Finally, we bound the centre–tail contribution
where
and
We can easily obtain the bounds
and
so again, by choosing *R* large enough we can ensure that
3.2

We therefore choose *R* such that conditions (3.1) and (3.2) are satisfied. Then we let be a sufficiently close analytic perturbation of *z*^{R}(*α*) which satisfies all the hypotheses of the lemma with *α*_{1}=*R* (so it also has vertical tangents at *α*=0,±*R*). It is not difficult to see that this can be done so that has the same sign as ∂_{α}*v*_{1}(*z*^{R}(*α*_{0})) for *α*_{0}=0,±*R*. Then will be the desired initial condition. ▪

### Lemma 3.2

*Let* *z*(*α*,*t*) *and* *w*(*α*,*t*) *be two analytic solutions of equation* (1.4) *with initial conditions* *z*_{0}(*α*) *and* *w*_{0}(*α*), *respectively. Let* *d*(*α*,*t*)=*z*(*α*,*t*)−*w*(*α*,*t*), *and for* let
*Assume* *and consider the energy*
*with* *S*(*t*)={*α*+*ic*:|*c*|<*ξ*(*t*)} *the strip of analyticity of the function* *d*(⋅,*t*). *Then there exists a polynomial* *P*(*x*,*y*,*q*,*r*,*s*) *such that*
*where the* *norm in the strip is defined as*

### Proof.

The proof appears in [4], Section 6, p. 940, eqn. (44). One only has to write the equation that *d*(*α*,*t*) satisfies, then apply the same estimates as in the local existence theorem from [4], and finally control the evolution of the norms of *z* and *w* in terms of the norms of *z*_{0} and *w*_{0} via their respective local existence theorems. ▪

We are now ready to prove our main result.

### Theorem 3.3

*There exists an analytic initial curve z(α,0) whose analytical extension is H*^{4} *on the boundary of its strip of analyticity, and also times −T< t*_{1}*<0<t*_{0}*<T and ε>0 such that the following hold. The corresponding solution of (1.4) exists for t∈(−T,T) in the class of analytic functions of α whose analytic extensions are H*^{4} *on the boundaries of their strips of analyticity, and it is a graph of a function of x for each t∈(t*_{1}*,t*_{0}*) (i.e. it is in the stable regime for these t) but not for t∈(t*_{1}*−ε,t*_{1}*)∪(t*_{0}*,t*_{0}*+ε) (i.e. it is in the unstable regime for these t). The solution develops vertical tangents at times t*_{1} *and t*_{0}.

### Proof.

For *δ*≥0, let *z*^{δ}(*α*) be a sufficiently close perturbation of with the same properties except that its tangents at *α*=0,±*R* are *δ* away from vertical (in the sense of for *α*_{0}=0,±*R*, while remains away from 0, uniformly in *δ*). It is not difficult to see that this can be chosen with a radius of analyticity away from 0, uniformly in *δ*. Then the solutions *z*^{δ}(*α*,*t*) of (1.4) corresponding to initial conditions *z*^{δ}(*α*) exist in the class of analytic functions on the interval (−*T*,*T*) for some *δ*-independent *T*>0.

Let *t*_{0}(*δ*) be the time of turnover of *z*^{δ}(*α*,*t*) near *α*_{0}=0, and let *t*_{1}(*δ*) be the time of turnover near *α*_{0}=±*R* (these exist if *z*^{δ}(*α*) is close enough to ). We have *t*_{0}(0)=*t*_{1}(0)=0, as well as *t*_{0}(*δ*)>0>*t*_{1}(*δ*) for sufficiently small *δ* and , due to lemma 3.2. Choosing *δ* such that −*t*_{1}(*δ*),*t*_{0}(*δ*)<*T* yields the desired initial condition *z*(*α*,0)=*z*^{δ}(*α*). ▪

## 4. Further discussion

We now discuss a couple of remarks related to theorem 3.3. We plan to address these in the future.

The first is the question of narrowing the gap between theorems 1.1 and 1.2, which was our initial motivation. Based on our numerical simulations, we propose the following conjecture.

### Conjecture 4.1

There exists a solution *f*(*x*,*t*) of equation (1.5) that has and turns over in finite time.

We present in figure 4 the supporting numerical evidence, with initial condition
If we reparametrize this curve as (*x*,*f*(*x*)), then .

A possible strategy to proving this conjecture, outlined in [17], is as follows. One can consider an approximate solution *w*(*α*,*t*), satisfying (with a small error err(*α*,*t*))
4.2
Next, subtracting equation (4.1) from (1.4) and defining *d*(*α*,*t*)=*z*(*α*,*t*)−*w*(*α*,*t*), one can write a system of equations for *d*(*α*,*t*) in such a way that only *d*(*α*,*t*) *w*(*α*,*t*) and *err*(*α*,*t*) appear, and *z*(*α*,*t*) does not. We recall that *w*(*α*,*t*) is explicit since it is the numerically calculated function and ∂_{α}*w*_{1}(*α*,*T*)<0 for some explicit (*α*,*T*). In a similar manner as in the local existence theorem from [4], a stability theorem follows, and this gives explicit bounds on the evolution of some norms of *d*(*α*,*t*), in particular it controls the norm of ∂_{α}*d*_{1}(*α*,*t*). The bounding of the constants of the stability theorem can be done either via traditional pencil-and-paper means, or using rigorous computer-assisted bounds and interval arithmetics, or a combination of both. Finally, if is small enough, then ∂_{α}*z*_{1}(*α*,*T*)=∂_{α}*w*_{1}(*α*,*T*)+∂_{α}*d*_{1}(*α*,*T*)<0. Note that, as opposed to lemma 3.2, we need to have good enough bounds on the constants that appear in the inequality, not just to know their existence. This makes the problem considerably harder.

### Remark 4.2

We believe that the constant 50 is not sharp and can be improved with further numerical search and better estimates.

### Remark 4.3

Showing existence of a stability shift in the other direction (stable unstable stable) seems harder, even though we have numerics that exhibit that behaviour. One has to produce an initial condition that first turns over and then recoils back, so in contrast with theorem 3.3, the solution lives in the unstable regime during the ‘middle’ time interval of its evolution.

A similar strategy as outlined above could work. After finding a numerical guess for a such solution, one can try to show via a stability theorem and computer-assisted estimates that close to this guess there exists a true solution of (1.4) which still exhibits this phenomenon.

## Competing interests

We declare we have no competing interests.

## Funding

D.C. and J.G.S. were partially supported by the grant no. MTM2011-26696 (Spain) and ICMAT Severo Ochoa project SEV-2011-0087. A.Z. acknowledges partial support by NSF grant nos. DMS-1056327 and DMS-1159133.

## Footnotes

One contribution of 15 to a theme issue ‘Free boundary problems and related topics’.

- Accepted February 9, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.