On an asymptotic model for free boundary Darcy flow in porous media

We provide a rigorous mathematical study of an asymptotic model describing Darcy flow with free boundary in a low amplitude/large wavelength approximation. In particular, we prove several well-posedness results in critical spaces. Furthermore, we also study how the solution decays towards the flat equilibrium.


Introduction
Since the pioneer works of Boussinesq [3], the derivation and study of asymptotic models for free boundary flows (usually, the water waves problem) are a hot research area (the interested reader can refer to [30]).
In this paper we study an asymptotic model for the intrusion of water into oil sand. This is known as the Muskat problem [31]. The study of the (full) Muskat problem has received a lot of attention in the last years, and, as a consequence, there is a large literature available (we refer to [26] for a recent survey of the available results). More precisely, this work considers the following one-dimensional equation   where Equation (AD ν ) was derived by the authors as an asymptotic model for the Darcy flow in porous media under the assumption that the amplitude over the wavelength (this ratio is known as the steepness) is small [27] (see also [9]).
The Bond number ν ≥ 0 represents a ratio between capillarity and gravitational effects. In the gravity driven case (when ν = 0) the equation (AD ν ) reads as   Although (AD ν ) and (AD 0 ) seem semilinear fourth and second order nonlocal PDE respectively (due to terms like f ∂ 4 x f and f ∂ 2 x f ), the commutator structure of the nonlinearity implies that they are a quasilinear third and a fully nonlinear first order nonlocal PDE respectively (see Lemma 3.1 below).
There are several motivations to study asymptotic models of free boundary Darcy flow. One of them arises from computational reasons. The idea is then, to simulate the asymptotic model to obtain an accurate description of the full problem at a lower computational cost. In that regards, let us briefly emphasize that the one-phase Muskat problem reveals itself as somehow harder (computationally speaking) than the two-phase Muskat problem. Let us try to briefly explain the reasons for this "paradoxical" fact. In the case when there are two fluids, the gravity driven Muskat problem reads as the following single nonlocal pde [13] ∂ t f = p.v.
Simplified models for this case where provided (following heuristic ideas) by Córdoba, Gancedo & Orive [15] (see also [29]). Remarkably, when the one phase Muskat problem is considered, the previous pde has to be modified and one is forced to study the following system of a nonlocal pde and an integral equation [12] ∂ t f (x) = p.v.
Another reason is the possibility of finding new finite time singularity scenarios. In particular, for the two-phase Muskat problem, Castro, Córdoba, Fefferman, Gancedo & López-Fernández [5] proved the existence of turning waves, i.e. interfaces that can be parametrized as a smooth graphs at time t = 0 but that become smooth curves that cannot be parametrized as graphs after a finite time (see also [2,[16][17][18]24]). We observe that these turning waves are interfaces such that there exists T 1 such that In the case of the one-phase Muskat problem, Castro, Córdoba, Fefferman, Gancedo [6] proved the existence of curves that self intersect in finite time in what is called a splash singularity (see also [19,[21][22][23]. It remains as an interesting open problem whether the Muskat problem can have a cusp singularity, i.e. a singularity where the slope and curvature of the interface blows up while the interface remains a graph (see [4,8,10,11,25,28] for global existence results). In that regards, we are optimistic enough to think that such an scenario should be simpler to prove (or discard) in an asymptotic model rather than the full problem. Taking this into consideration, the results of this paper then give conditions that excludes finite time blow ups of cusp and also turning type.

Functional spaces
The space domain considered in the present article is the one-dimensional torus. i.e. T = R 2π Z. The domain T can also be understood as the interval [−π, π] endowed with periodic boundary conditions. Let f (x) denote a L 2 function on T. Then, its Fourier transform is an ℓ 2 (Z) sequence defined aŝ for any n ∈ Z, with inverse Fourier transform Then, we define the L 2 −based (homogeneous) Sobolev spacesḢ α (T) In the present work we will need to define Sobolev spaces with fractional derivatives and integrability indexes different from two, i.e. spacesẆ s,p (T) for s ∈ R and p ∈ [1, ∞]. We provide here a characterization of such spaces using using the theory of Littlewood-Paley (see Appendix A). If we consider the dyadic block △ q (we refer again to Appendix A for a precise definition) we can hence define the semi-norm 1 in the sense that there exists a positive absolute constant C > 0 s.t. 1 C u Ḣ s (T) u Ẇ s,2 (T) C u Ḣ s (T) . Let us recall the following embedding, valid for mean-free distributions on the one-dimensional torus: and in particular the following inequality holds true for any u ∈Ḣ 1/2 ∩ L p , Similarly, we define the (homogeneous) Wiener spacesȦ α (T) and Wiener spaces with weight aṡ respectively. Obviously,Ȧ α =Ȧ α 0 and the functions in A s ν are analytic in the complex strip We observe that the equation (AD ν ) conserves the average of a solution f , whence w.l.o.g. we can assume this average to be zero. Thus, from this point onwards we identify the spacesḢ s (T) and H s (T),Ȧ s (T) and A s (T) andȦ s ν (T) and A s ν (T).

A word on scaling
We observe that, when ν = 0, equation (AD ν ) is left invariant by the scaling Remarkably, this scaling is the same as in the full Muskat problem (see [20]). Then, we note that the following spaces are critical for this scaling Let us now consider the equation (AD ν ) for ν > 0. In such setting (AD ν ) reads, when expanded, as There exists indeed no scale invariance satisfied by the equation (1.3), despite this fact we can rewrite (1.3) as This result is a sort of Cauchy-Kovalevsky Theorem. However, we emphasize that the solution is less regular in time (merely L ∞ instead of continuous) but maintains its original strip of analyticity.
Furthermore, we are able to prove a result establishing the decay of certain norms of the solution: This result seems the analog of the maximum principles known for the full Muskat problem [14]. In particular, when compared to the results in [14,18], it seems that the model equation (AD 0 ) is less stable than the full Muskat problem.
Once the local solution for analytic data is known, we turn our attention to the global solution for initial data satisfying certain size restrictions in critical spaces. Then we state our well-posedness result for Wiener class initial data: 3. Let f 0 ∈ A 1 be the initial data for (AD 0 ). Assume that then the Cauchy problem (AD 0 ) is globally well-posed and admits a unique solution We observe that this is a global existence and decay for initial data in a critical space.
The next result we prove in the present manuscript is a global existence result for initial data in H  [20] that it possible to construct global solutions for the 2D (full) Muskat problem when the initial data is small in The following theorem establishes a similar result for the evolution equation (AD 0 ); the Cauchy problem (AD 0 ) is globally well-posed and admits a unique solution f such that for any T > 0 The methodology used in order to prove Theorem 2.4 differs completely from the one used in [20], since (AD 0 ) presents a nonlinearity of polynomial type we will be able to use tools characteristic of the paradifferentail calculus, we refer the reader to [1, Chapter 2]; we will hence decompose the nonlinearity in an infinite sum of elementary packets on which we will be able to highlight some nontrivial regularizing properties of the equation (AD 0 ).

Results for the gravity-capillary driven case (AD ν )
We are able to extend Theorem 2.3 to the case with surface tension: Theorem 2.5. Let ν > 0 be a fixed parameter and f 0 ∈ A 1 be the initial data for (AD ν ). Assume that then the Cauchy problem (AD ν ) is globally well-posed and admits a unique solution Remark 2.6. We observe that the size condition is ν−independent.
The following result is a result analogous of the one stated in Theorem 2.4 for the system (AD ν ) when ν > 0:

Discussion
In this paper we prove several well-posedness results for an asymptotic model of free boundary Darcy flow. Most of them are global existence results in scale invariant spaces. In that regards, our results should be understood as nonlinear stability results rather than linear stability (even if some size conditions are imposed on the initial data).
These results excludes the existence of turning waves singularities but leave open the door to cusp singularities. The occurrence of such behavior will be the object of future research.

Proof of Theorem 2.1
We look for a solution of the form where λ = λ( f 0 ) will be chosen below. Then, the existence of solution is reduced to the summability of a series where each term satisfy a linear problem. A similar idea can be tracked back to the works of Oseen [32] (see also [9]). Indeed, matching the appropriate terms, we find that f (ℓ) satisfies with initial data Then and we can solve recursively for the other f (ℓ) . In particular, using that the solution of is given byû we find that Thus, using Thus, we obtain that We now fix 1 ≪ R and consider the partial sum .
We estimate Applying Tonelli's theorem, we find that and Thus, We define the numbers These numbers satisfy Then, we prove by induction that where C ℓ are the Catalan numbers. It is known that the Catalan numbers grow like 4 ℓ . We also observe that As a consequence, we have the bound Thus, we have a bound This bound is independent of R and then, using Weierstrass M-theorem in the space L ∞ (0, T ; A 1 1 ), we can pass to the limit as R → ∞ and obtain Then, using the Cauchy theorem for the product of series, the limit is a mild solution of our problem

Proof of Proposition 2.2
We start this section with an auxiliary Lemma that establishes an integral formula for the nonlinear term Proof of Lemma 3.1. First, we observe that Thus, we have that where tan(z/2) dzdy, tan(z/2) dzdy, We have that Using that, for zero mean functions, together with (3.3), we find that Thus, collecting every term, we conclude the result.

Proof of Proposition 2.2.
Define X t such that We observe that M(t ) is a.e. differentiable. To see that take two t ≤ s and assume that M(t ) ≥ M(s), then As a consequence, the boundedness of ∂ t f L ∞ is a sufficient condition for M(t ) to be Lipschitz. Now Rademacher Theorem gives us the a.e. differentiability of M. Furthermore, one can prove that the derivative verifies the following almost everywhere in time. Using Lemma 3.1, we find that If we define x t such that we can repeat the previous step and conclude the result.

Proof of Theorem 2.3
Step 1: a priori estimates We start the proof providing the appropriate a priori estimates in A 1 . We have that In Fourier variables, we have that We observe that p = 0 if and only if 0 < |k| < |m|. We find that |p(k, m)| |k|2|k||k − m| ≤ 2|m| 2 |k − m|. Thus, As a consequence, we find the inequality Step 2: Approximated solutions and passing to the limit To construct solutions we consider the regularized problem where the initial data is localized in Fourier space. In other words, given f 0 , we consider (AD 0 ) with the initial data f N 0 where Then, invoking Theorem 2.1, there exists a local solution. For this local solution we can apply the previous energy estimates and, consequently, we can pass to the limit.

Step 3: Uniqueness
We argue by contradiction: let's assume that there exist two different solutions f 1 and f 2 starting from the same initial data. Then, one can prove that the difference g = f 1 − f 2 satisfies d dt Then, one concludes the uniqueness using Gronwall's inequality.

Step 4: Decay
We have that d dt Recalling (3.1), we find that As a consequence, we find the inequality Using the Poincaré-like inequality we find the exponential decay

Proof of Theorem 2.4
Step 1: a priori estimates We consider f to be a (space-time) smooth solution of the periodic problem (AD 0 ). For such solutions we want to prove the following global energy estimate then the corresponding solution f of (AD 0 ) satisfies: In order to outline the idea behind the proof of Proposition 3.2 we shall, at first, try to estimate globally the L 2 (T) norm of a suitable solution. In that regards state and prove the following lemma: Proof of Proposition 3.3. Let us multiply equation (AD 0 ) for f and integrate by parts. Then, we deduce where N g 1 , g 2 , g 3 = − T g 1 (H − 1) ∂ x g 2 (H + 1) ∂ x g 3 dx. (3.8) We rely now on the following technical lemma, which entails the commutation property mentioned above and that concludes the proof of Proposition 3.3. We divide the proof of Proposition 3.2 in several steps: Step 1.1: dyadic estimates.
Let us apply the truncation operator △ q to the equation (AD 0 ), multiply the resulting equation for △ q f and integrating in space we obtain: Using now Bony decomposition (A.2) we can say We have to estimate these eight terms. The terms T A 1 and T B 1 are the less regular. We will use the commutation properties proved in Lemma 3.4 to estimate such terms. First of all we remark that where the trilinear operator N is defined in (3.8). We can hence use Lemma 3.4 and (A.3) in order to deduce for any s 3/2.
Next we handle the remainder terms R A q , R B q , indeed Since H 1 2 (T) → L p (T) for any p ∈ [2, ∞) we can say that Moreover due to the localization in the Fourier space and the fact that q ′ > q − 4 we deduce (3.11) We can deduce as well the following estimate for the term R A for any 0 < δ ≪ 1.
A very similar procedure allows us to deduce the following bound We turn now our attention to the element T A 2 . Using Lemma A.1 we obtain Again with simple computations similar to the ones performed above we can conclude the following bound: With an analog procedure we deduce the bounds The estimates can be easily obtained since the terms composing the elements T A 3,q , T B 3,q are all localized in dyadic annuli, and hence such terms are more regular than the ones estimated above. We note that where we have used (3.10), (3.11), (3.13), (3.15), (3.17) and (3.19). Now, multiplying the above equation for 2 2qs and summing in q ∈ Z, we obtain (3.20) Using (3.10), (3.12), (3.14), (3.16), (3.18) and (3.19) instead we deduce Step 1.2: H 3/2 estimates We can now simply consider the estimate (3.20) in which we set s = 3/2 and in order to deduce 1 2 d dt f Step 1.5: proof of (3.5) We multiply AD 0 for ∂ t Λ 2 f and integrate in T × (0, T ) obtaining the following integral equality (here L p (3.25) Using the identity Λ 2 = −∂ 2 x and integrating by parts we can easily deduce the following identity = I 1,1 + I 1,2 , We hence obtained that The term C 1 can be estimated using Lemma 3.4 and we deduce The term C 2 can be estimated performing computations very similar to the ones used in order to prove Lemma 3.4, in particular we obtain for any ǫ > 0 What remains hence to estimate is the term I 1,1 + I 2,1 , we will perform the estimates for the term I 1,1 only being the ones for I 2,1 very similar (here 2 < p ε = ε −1 < ∞) (3.28) We hence plug the estimates (3.26), (3.27)(3.28) in (3.25) obtaining We fix δ = 1/4 and we use the main estimate of Proposition 3.2, i.e. estimate (3.4) in order to deduce the bound 4 Gravity-capillarity driven system (AD ν )

Proof of Theorem 2.5
Step 1: a priori estimates We only prove the a priori estimates, since the approximation procedure is standard. We have that As before, In Fourier variables, we have that Again p = 0 if and only if 0 < |k| < |m|. We find that Thus, using the interpolation inequality (valid for 0 ≤ r ≤ 3) we find that As a consequence, we find the inequality and we conclude the result.
Step 2: Decay From here we conclude using a Poincaré type inequality.

Proof of Theorem 2.7
Step 1: a priori estimates We start with a Proposition giving conditions that guarantee the existence of an L 2 energy balance: Let us consider f a smooth solution of (AD ν ) in T, and let suppose that, for any t ∈ [0, T ] then the following inequality holds true for any t ∈ [0, T ], Proof. As it was done in section 3, we simply perform L 2 energy estimates for (AD ν ) and we study the particular commutation properties of the nonlinearity of (AD ν ). Let us recall that whence multiplying (AD ν ) for f and integrating in T we deduce the following where N is defined in (3.8) and M is defined as We rely now on the following lemma, which can be understood as the analog of Lemma we deduce that which closes the estimates.
We want now to prove global H 2 estimates for (AD ν ) stemming from small initial data. The result we prove is the following one: then for any t 0 the following estimate holds true Proof. Again as in the proof of Proposition 3.2 we use the tools of paradifferential calculus illustrated in Appendix A. We apply the dyadic truncation △ q to the left of (AD ν ), we multiply the resulting equation for △ q f and integrate in x We can use the estimates performed in Proposition 3.2 (see in particular the r.h.s. of (3.21)) in order to deduce We can hence now focus on the purely nonlinear part which is characteristic of (AD ν ) when ν > 0, i.e. −ν C q + D q . As it was done in the proof of Proposition 3.2 we can use Bony decomposition A.2 in order to decompose C q and D q as We remark now that Whence we can apply Lemma 4.2 and (A.3) in order to deduce the bound Next we focus on the remainder terms, as in the proof of Proposition 3.2 we can handle R C q and R D q hence we provide an explicit computation for the former only. Using Hölder inequality, Sobolev embeddings, (A.3) and interpolation of Sobolev spaces it is possible to deduce the following estimate In the above estimates indeed which is ℓ 2 being c q q ∈ ℓ 1 . Similar computations holds for the term R D q , from where we deduce that Next we study the term T C 2,q . Using Lemma A.1 we obtain △ q ′ Λ 3/2 f L 2 and as a consequence of 2 q ′ ∼ 2 q for q − q ′ 4 we find that A similar estimate holds for T D 2,q . The terms T C 3,q , T D 3,q enjoy a analog bounds being composed by linear combinations of functions localized in annuli. We hence deduced that (4.6)

A Elements of Littlewood-Paley theory.
A tool that will be widely used all along the paper is the theory of Littlewood-Paley, which consists in doing a dyadic cut-off of the frequencies. Let us define the (homogeneous) truncation operators as follows: where u ∈ D ′ T 3 andû n are the Fourier coefficients of u. The function ϕ is a smooth function with compact support such that and such that for all t ∈ R, q∈Z ϕ 2 −q t = 1.
Let us define further the low frequencies cut-off operator

A.1 Paradifferential calculus.
The dyadic decomposition turns out to be very useful also when it comes to study the product between two distributions. We can in fact, at least formally, write for two distributions u and v Paradifferential calculus is a mathematical tool for splitting the above sum in three parts The following almost orthogonality properties hold and hence we will often use the following relation In the paper [7] J.-Y. Chemin and N. Lerner introduced the following decomposition which will be very useful in our context △ q (uv) = S q−1 u △ q v + |q−q ′ | 4 where the commutator △ q , a b is defined as Let us define the following distribution H (x) = F −1 ĥ n n (x) .
A sufficient condition for H to be well defined is that ĥ n n ∈ ℓ 2 , if this is the case H will result to be a L 2 function.