Analytical behavior of two-dimensional incompressible flow in porous media

ANALYTICAL BEHAVIOR OF 2-D INCOMPRESSIBLE FLOW IN POROUS MEDIA DIEGO C´ ORDOBA, FRANCISCO GANCEDO AND RAFAEL ORIVE Abstract. In this paper we study the analytic structure of a two-dimensional mass balance equation of an incompressible fluid in a porous medium given by Darcy’s law. We obtain local existence and uniqueness by the particle-trajectory method and we present different global existence criterions. These analytical results with numerical simulations are used to indicate nonformation of singularities. Moreover, blow-up results are shown in a class of solutions with infinite energy. 1. Introduction The dynamics of a fluid through a porous medium is a complex and not thoroughly understood phenomenon [3, 16]. The purpose of the paper is to study a nonlinear two-dimensional mass balance equation in porous media and the conditions of formation of singularities using analytical results and numerical calculations. In real applications, one might be interested in the transport of a dissolved contaminant in porous media where the contaminant is convected with the subsurface water [17]. For example, one is usually interested in the time taken by the pollutant to reach the water table below. Such flows also occur in artificial recharge wells where water and (or) chemicals from the surface are transported into aquifers. In this case, the chemicals may be the nutrients required for degradation of harmful polluting hydrocarbons resident in the aquifer after a spillage. We use Darcy’s law to model the flow velocities, yielding the following relationship between the liquid discharge (flux per unit area) v= (v1, v2) and the pressure gradient v=−k(∇p+gγρ), where kis the matrix position-independent medium permeabilities in the different directions respectively divided by the viscosity, ρis the liquid density, gis the acceleration due to gravity and the vector γ= (0,1). While the Navier-Stokes equation and the Stokes equation are both microscopic equations, Darcy’s law gives the macroscopic description of a flow in a porous medium [3]. The free boundary problem given by an incompressible 2-D fluid through porous media with two different constant densities and viscosities at each side of the interface is studied in [18, 1] (see references therein). Here, we analyze the dynamics of the density function ρ=ρ(x1, x2, t) with a regular initial data ρ0=ρ(x1, x2,0). The mass balance equation is given by φDρ Dt =φ∂ρ ∂t +v· ∇ρ= 0, where φdenotes the porosity of the medium. To simplify the notation, we consider k=g=φ= 1. Thus, our system of a two-dimensional mass balance equation in porous media (2DPM) is written as Dρ Dt = 0,(1.1) v=−(∇p+γρ).(1.2) We close the system assuming incompressibility, i.e., (1.3) divv= 0, Date: May 31, 2006. Key words and phrases. Flows in porous media 2000 Mathematics Subject Classification. 76S05, 76B03, 65N06. The authors was partially supported by the grant MTM2005-05980 of the MEC (Spain) and S-0505/ESP/0158 of the CAM (Spain). The third author was partially supported by the grant MTM2005-00714 of the MEC (Spain). 1 2 D. C ´ ORDOBA, F. GANCEDO AND R. ORIVE therefore there exists a stream function ψ(x, t) such that (1.4) v=∇⊥ψ≡−∂ψ ∂x2 ,∂ψ ∂x1. Computing the curl of equation (1.2), we get the Poisson equation for ψ (1.5) −∆ψ=∂ρ ∂x1 . A solution of this equation is given by the convolution of the Newtonian potential with the function ∂x1ρ (1.6) ψ(x, t) = −1 2πZR2 ln|x−y|∂ρ ∂y1 (y, t)dy, x ∈R2. Thus, the velocity vcan be recovered from ψby the operator ∇⊥by the two equivalent formulas v(x, t) = ZR2 K(x−y)∇⊥ρ(y, t)dy, x ∈R2,(1.7) v(x, t) = PV ZR2 H(x−y)ρ(y, t)dy −1 2(0, ρ(x)) , x ∈R2,(1.8) where the kernels K(·) and H(·) are defined by K(x) = −1 2π x1 |x|2and H(x) = 1 2π−2x1x2 |x|4,x2 1−x2 2 |x|4.(1.9) Differentiating the equation (1.1), we obtain the evolution equation for ∇⊥ρ≡−∂ρ ∂x2 ,∂ρ ∂x1 (1.10) which is given by (1.11) D∇⊥ρ Dt = (∇v)∇⊥ρ. Taking the divergence of the equation (1.2), we get −∆p=∂ρ ∂x2 and the pressure can be obtained as in (1.6). The objective of this work is to analyze the behavior of the solutions of the system 2DPM (1.1)- (1.3). First, we present the existence of singularities in a class of solutions with infinite energy in 2DPM (see Proposition 2.2 and Remark 2.3). In the case of solutions with regular initial data and finite energy, we get local well-posedness using the classical particle trajectories method. We illustrate a criterion of global existence solutions via the norm of the bounded mean oscillation space of (1.10). A similar result is known in the three dimensional Euler equation (3D Euler) [2]. Also, using the geometric structure of the level sets of the density (where ρis constant) and the nonlinear evolution equations of the gradient of the arc length of the level sets, we establish that no singularities are possible under not very restrictive conditions. This result is comparable to the 3D Euler equations [7] and to the two-dimensional quasi geostrophic equation (2DQG) [8]. Applying these criterions, we find no evidence of formation of singularities in our numerical simulations. The paper is organized as follows. In Section 2 we study the analytical behavior of solutions with infinite energy. In Section 3 we prove the existence and uniqueness for the 2DPM, show a characterization of formation of singular solutions and we present geometric constraints on singular solutions. Finally, in Section 4 we illustrate two numerical examples in which the analytical results are applied to show nonsingular solutions. ANALYTICAL BEHAVIOR OF 2DPM 3 2. Singularities with infinite energy Let the stream function ψbe defined by (2.1) ψ(x1, x2, t) = x2f(x1, t) + g(x1, t). Note that under this hypothesis the solution of (1.1)–(1.3) has infinite energy. We reduce the equations (1.1)–(1.3) to other system with respect to the functions fand g. From (1.5) the density, apart from a constant, satisfies (2.2) ρ(t, x1, x2) = −x2 ∂f ∂x1 (x1, t)−∂g ∂x1 (x1, t) = −x2fx1−gx1, and, by (1.4), vverifies (2.3) v(t, x1, x2) = −f(x1, t), x2 ∂f ∂x1 (x1, t) + ∂g ∂x1 (x1, t)= (−f, x2fx1+gx1). Therefore, the system (1.1)–(1.3) under the hypothesis (2.1) is equivalent to (fx)t=ffxx −(fx)2,(2.4) (gx)t=fgxx −fxgx.(2.5) (Here and in the sequel of the section, we denote with subscript the derivatives with respect to x.) We note the non-linear character of the first equation. Thus, our study of formation of singularities is concentrated in the solutions of (2.4). The function gdepends implicitly on fin equation (2.5). Now, we show that the system (2.4) and (2.5) is local well posed in the Sobolev spaces Hk 0(0,1). Lemma 2.1. Let f0=f(x, 0) and g0=g(x, 0) satisfy f0 x, g0 x∈Hk 0(0,1) with k≥1. Then, there exists T > 0such that fx, gx∈C1([0, T]; Hk 0(0,1)) are the unique solution of (2.4)–(2.5). Proof. By (2.4) and integrating by parts, we have 1 2 d dtkfxk2 L2=Z1 0 fxffxx −Z1 0 f3 x=−3 2Z1 0 f3 x≤CkfxkL∞kfxk2 L2≤Ckfxk3 H1 0. Analogously, 1 2 d dtkfxxk2 L2=−Z1 0 f2 xxfx−Z1 0 fxxffxxx =−1 2Z1 0 f2 xxfx≤CkfxkL∞kfxxk2 L2≤Ckfxk3 H1 0. We can repeat for all k≥1 and we obtain 1 2 d dtkfxk2 Hk 0≤Ckfxk3 Hk 0. Integrating in time, we get kfxkHk 0≤kf0 xkHk 0 1−Ctkf0 xkHk 0 . On the other hand, by (2.4) and integrating by parts, we have for gxthe following inequalities 1 2 d dtkgxk2 L2=Z1 0 gxgxxf−Z1 0 g2 xfx=−3 2Z1 0 g2 xfx≤ kfxkL∞kgxk2 H1 0 and 1 2 d dtkgxxk2 L2=Z1 0 gxxgxxxf−Z1 0 gxxgxfxx =−1 2Z1 0 g2 xxfx−Z1 0 gxxgxfxx ≤ kfxkL∞kgxxk2 L2+kgxxkL2kgxkL∞kfxxkL2≤ kfxkH1 0kgxk2 H1 0. Thus, we obtain using Gronwall’s Lemma kgxk2 Hk 0≤ kg0 xk2 Hk 0exp CZt 0 kfxkHk 0ds and we have existence up to a time T=T(kf0 xkHk 0). 4 D. C ´ ORDOBA, F. GANCEDO AND R. ORIVE In order to prove the uniqueness, let fx(x, t) = hx(x, t)−kx(x, t), with hx, kxtwo solutions of (2.4) with the same initial data f0 x. Since hx, kxsatisfy (2.4) and integrating by parts, we have 1 2 d dtkfxk2 L2=Z1 0 fx(hhxx −kkxx)−Z1 0 fx(h2 x−k2 x) =Z1 0 fxhfxx +Z1 0 fxfkxx −Z1 0 (fx)2(hx+kx) =−1 2Z1 0 (fx)2hx+Z1 0 fxfkxx −Z1 0 (fx)2(hx+kx). Thus, we get 1 2 d dtkfxk2 L2≤ kkxxkL2kfxkL2kfkL∞+C(khxkL∞+kkxkL∞)kfxk2 L2 ≤C(khxkH1 0+kkxkH1 0)kfxk2 L2 and using Gronwall’s Lemma it follows hx=kx. Finally, we conclude the uniqueness of gxsince (2.5) is a linear differential equation . The following result shows that the solution of (2.4) blows up in finite time under certain conditions on the initial data. Proposition 2.2. Let fxbe a solution of (2.4) with initial data satisfies f0 x∈H2 0(0,1) and minxf0 x< 0. Then, kfxkL∞blows up in finite time T=−1/minxf0 x. Proof. By the local existence result, we have fx∈C1([0, T ]; H2)⊂C1([0, T]×[0,1]). We consider the application m: [0, T]→Rdefined by m(t) = minxfx(x, t) = fx(xt, t). By Rademacher Theorem, it follows that mis differentiable at almost every point. First, we calculate the derivative of mas in [5, 9]. Let sbe a point of differentiability of m(t), then for τ > 0 m0(s) = lim τ→0 m(s+τ)−m(s) τ= lim τ→0 fx(xs+τ, s +τ)−fx(xs, s) τ = lim τ→0 fx(xs+τ, s +τ)−fx(xs, s +τ) τ+fx(xs, s +τ)−fx(xs, s) τ. Since fx(x, s +τ) reaches the minimum at the point xs+τ, we obtain m0(s)≤lim τ→0 fx(xs, s +τ)−fx(xs, s) τ=fxt(xs, s). We compute the derivative with a sequence of negative τ < 0 and, by the sign of τ, we get the opposite inequality and we conclude that m0(s) = fxt(xs, s) almost everywhere. We replace xfor xsin (2.4) and yields m0(s) = −f2 x(xs, s) = −(m(s))2, due to fxx(xs, s) = 0, and the proof follows. Remark 2.3. There are other blow-up results with an initial data of lower regularity. In particular, we consider f0 x∈H1 0and assuming that Z1 0 f0 x≤0. Thus, by (2.4), we have d dt Z1 0 fx=Z1 0 ffxx −Z1 0 (fx)2=−2Z1 0 (fx)2≥ −2Z1 0 fx2 . Defining c(t) = Zfx, ANALYTICAL BEHAVIOR OF 2DPM 5 and integrating, we get c(t)≤c(0) 1+2tc(0). Then, c(t)blows up for c(0) <0. In the case c(0) = 0, we have c0(t)<0for all t > 0, therefore, c(t)also blows up. Remark 2.4. Let x1=xtbe the point such that fx(xt, t) = min xfx(x, t), and consider x2= 1 −gx(xt, t) fx(xt, t). Then, by (2.2), ρ(x1, x2, t) = −fx(xt, t)blows up in finite time by Proposition 2.2. Analogously, v defined in (2.3) blows up in finite time. 3. Analysis of 2DPM with finite energy 3.1. Local existence of 2DPM. We derive a reformulation of the system as an integro-differential equation for the particle trajectories. Given a smooth field v(x, t), the particle trajectories Φ(α, t) satisfy dΦ dt (α, t) = v(Φ(α, t), t),Φ(α, t)|t=0 =α.(3.1) The time-dependent map Φ(·, t) connects the Lagrangian reference frame (with the variable α) to the Eulerian reference frame (with the variable x). It is well known (Section 2.5 in [15]) that the equation (1.11) implies the following formula ∇⊥ρ(Φ(α, t), t) = ∇αΦ(α, t)∇⊥ρ0(α), where ∇⊥ρ0is the orthogonal gradient of the initial density. This last equality shows us that the orthogonal gradient of the density is stretched by ∇αΦ(α, t) along particle trajectories. We rewrite (1.7) as (3.2) v(Φ(α, t), t) = ZR2 K(Φ(α, t)−Φ(β, t)) ∇αΦ(β, t)∇⊥ρ0(β)dβ. We consider (3.1) as an ODE on a Banach space and using Picard Theorem the local in time existence follows. This is proved analogously like the existence and uniqueness of solutions to the inviscid Euler equation (see Section 4.1 in [15]). In fact, we consider ∇⊥ρ0∈Cδ(R2), δ∈(0,1). Let Bbe the Banach space defined by B=Φ : R2→R2such that |Φ(0)|+|∇αΦ|0+|∇αΦ|δ<∞, where |·|0is the L∞-norm and |·|δis the H¨older semi-norm. Define OM, the open set of B, as OM=Φ∈B|inf α∈R2det ∇αΦ(α)>1 2and |Φ(0)|+|∇αΦ|0+|∇αΦ|δ< M. The mapping v(Φ), defined by (3.2), satisfies the assumptions of the Picard theorem, i.e., vis bounded and locally Lipschitz continuous on OM. As a consequence, for any M > 0 there exists T(M)>0 and a unique solution Φ∈C1((−T(M), T(M)); OM) to the particle trajectories (3.1, 3.2). Remark 3.1. The 2DPM (1.1)–(1.3) has quantities conserved in time, the Lpnorm of ρfor 1≤p≤ ∞, i.e., (3.3) kρ(t)kp=kρ0kp,∀t > 0,1≤p≤ ∞. The velocity is obtained from ρby (1.8). These operators are singular integrals with Calder´on-Zygmund kernels (see [19]). Then for 1< p < ∞the Lpnorm of the velocity is bounded for any time t > 0. 6 D. C ´ ORDOBA, F. GANCEDO AND R. ORIVE 3.2. Blow-up criterion. In order to estimate the growth of the Sobolev norms we use the space of functions of bounded mean oscillation (BMO) (see Chapter IV in [19] for an introduction of this function space). Theorem 3.2. Let ρbe the solution of equation (1.1)–(1.3) with initial data ρ0∈Hs(R2)with s > 2. Then, the following are equivalent: (A) The interval [0,∞)is the maximal interval of Hsexistence for ρ. (B) The quantity (3.4) ZT 0 k∇ρkBMO(t)dt < ∞ ∀T > 0. Proof. We denote the operator Λsby Λs≡(−∆)s/2. Since the fluid is incompressible, we have for s > 2 1 2 d dtkΛsρk2 L2=−ZR2 ΛsρΛs(v∇ρ)dx =−ZR2 Λsρ(Λs(v∇ρ)−vΛs(∇ρ)) dx, ≤CkΛsρkL2kΛs(v∇ρ)−vΛs(∇ρ)kL2. Using the following estimate (see [13]) kΛs(fg)−fΛs(g)kLp≤ck∇fkL∞kΛs−1gkLp+kΛsfkLpkgkL∞1< p < ∞, we obtain for p= 2 (3.5) 1 2 d dtkΛsρk2 L2≤C(k∇vkL∞+k∇ρkL∞)kΛsρk2 L2. Integrating, we get for any t≤T (3.6) kΛsρkL2≤ kΛsρ0kL2exp CZT 0 (k∇vkL∞+k∇ρkL∞)!. Now, we use the following inequality given in [14]: Let f∈Ws,p with 1 <p<∞and s > 2/p, then, there exists a constant C=C(p,s) such that (3.7) kfkL∞≤C(1 + kfkBMO(1 + ln+kfkWs,p )), where ln+(x) = max(0,ln(x)). Therefore, for s > 2 we have k∇ρkL∞≤C(1 + k∇ρkBMO(1 + ln+k∇ρkHs−1)), and from (3.6), we obtain that (3.8) k∇ρkL∞≤C(1 + ln+(kρ0kHs)k∇ρkBMO ZT 0 (k∇vkL∞+k∇ρkL∞)dt). On the other hand, applying (3.7) for v∈Hs(R2), we have k∇vkL∞≤C(1 + k∇vkBMO(1 + ln+k∇vkHs−1)). Since vsatisfies (1.8) and the singular integrals are bounded operators in BMO (see [19]), we get k∇vkL∞≤C(1 + k∇ρkBMO(1 + ln+k∇ρkHs−1)), and, using (3.6), we obtain (3.9) k∇vkL∞≤C(1 + ln+(kρ0kHs)k∇ρkBMO ZT 0 (k∇vkL∞+k∇ρkL∞)dt). From (3.8) and (3.9), follows k∇vkL∞+k∇ρkL∞≤C(1 + ln+(kρ0kHs)k∇ρkBMO ZT 0 (k∇vkL∞+k∇ρkL∞)dt). Applying Gronwall’s inequality, we have ZT 0 (k∇vkL∞+k∇ρkL∞)dt ≤CT exp ln+(kρ0kHs)ZT 0 k∇ρkBMO dt!, and so (A) is a consequence of (B). ANALYTICAL BEHAVIOR OF 2DPM 7 Finally, due to the inequality k∇ρkBMO ≤ k∇ρkH1, we conclude that (A) implies (B). Remark 3.3. Using that k∇ρkBMO ≤Ck∇ρkL∞, we get a blow-up characterization for numerical simulations. 3.3. Geometric constraints on singular solutions. From equation (1.1) it follows that the level sets, ρ= constant, move with the fluid flow. Then ∇⊥ρ, defined in (1.10), is tangent to these level sets. For the 2DPM, the infinitesimal length of a level set for ρis given by |∇⊥ρ|and from (1.11), the evolution equation for the infinitesimal arc length is given by (3.10) D|∇⊥ρ| Dt =L |∇⊥ρ|. The factor L(x, t) is defined through by (3.11) L(x, t) = Dη·η, η 6= 0, 0, η = 0. where the direction of ∇⊥ρis denoted by (3.12) η=∇⊥ρ |∇⊥ρ| and D(x, t) is the symmetric part of the deformation matrix defined by (3.13) D= (Dij) = 1 2∂vi ∂xj +∂vj ∂xi. Now, we show a singularity criterion of 2DPM using the geometric structure of the level sets and mild hypotheses of the solutions. The theorem states below is analogous to 3D Euler [7] and to 2DQG [8]. Recall that ηis the direction field tangent to the level sets of ρdefined (3.12). Analogously to [8], a set Ω is smoothly directed if there exists δ > 0 such that (3.14) sup x∈ΩZT 0 k∇η(·, t)k2 L∞(Bδ(Φ(x,t)))dt < ∞, where Bδ(x) = {y∈R2:|x−y|< δ},Ω = {x∈Ω; |∇ρ0(x)| 6= 0}, and Φ is the particle trajectories map. We define Ω(t) = Φ(Ω, t) and OT(Ω) the semi-orbit, i.e., OT(Ω) = [ 0≤t≤T {t} × Ω(t). Theorem 3.4. If Ωis smoothly directed and (3.15) ZT 0 kRjρkL∞(t)dt < ∞, j = 1,2,∀T > 0, where Rjdenotes the Riesz transform in the direction xj, then sup OT(Ω) |∇ρ(x, t)|<∞. Remark 3.5. Using the Remark 3.3, the previous theorem illustrates that finite-time singularities are impossible in smoothly directed sets. Proof. We show a similar formula of the level–set stretching factor Ldefined in (3.11). We start by computing the full gradient of the velocity v. From formula (1.7) v(x) = ZR2 K(y)∇⊥ρ(x−y)dy, 8 D. C ´ ORDOBA, F. GANCEDO AND R. ORIVE we have ∇v(x) = ZR2 K(y) (∇y∇⊥ yρ)(x−y)dy. Take the integral as a limit as ε→0 of integrals on |y|> ε and integrate by parts. In this fashion, we obtain the formula (3.16) ∇v(x) = 1 2πPV ZR2∇yρ(x−y)⊗eydy |y|2−1 2   0 0 ∂ρ ∂x1 (x)∂ρ ∂x2 (x)  , where eyis the unit vector defined by ey=−2y1y2 |y|2,y2 1−y2 2 |y|2. By definition of ηin (3.12), we have η· ∇ρ= 0. Thus, computing we get (3.17) L(x) = 1 2πPV ZR2ey·η(x)η(x−y)·η⊥(x)}|∇⊥ρ(x−y)|dy |y|2. Let be φ∈C∞ c(R), φ≥0, supp(φ) include in [−1,1] and φ(s) = 1 in s∈[−1/2,1/2]. Consider r > 0 and decompose L(x) = I1+I2, with |I1| ≤ 1 2πZR2 (1 −φ(|y|2/r2))ey·η(x)∇⊥ρ(x−y)·η⊥(x)dy |y|2. Integrating by parts and using Cauchy-Schwartz inequality, we get I1≤C r2kρkL2≤C r2kρ0kL2. We have for any |y|< r |η(x−y)·η⊥(x)| ≤ |y|k∇ηkL∞(Br(x)). Applying this in the integral I2, we get |I2| ≤ Z|∇⊥ρ(x−y)|φ(|y|2/r2)dy |y|k∇ηkL∞(Br(x)). We integrate by parts and decompose Z|∇⊥ρ(x−y)|φ(|y|2/r2)dy |y|=Zρ(x−y)∇⊥η(x−y)φ(|y|2/r2)1 |y|dy =J1+J2+J3, where J1=Zρ(x−y)η(x−y)∇⊥(φ(|y|2/r2))dy |y|, J2=Zρ(x−y)∇⊥(η(x−y))φ(|y|2/r2)dy |y|, J3=Zρ(x−y)η(x−y)φ(|y|2/r2)(−y2, y1) |y|3dy, obtaining the following estimates |J1| ≤ ckρ0kL∞and |J2| ≤ crkρ0kL∞k∇ηkL∞(Br(x)). The J3term can be bounded using the identity J3=η(x)(−R2(ρ)(x), R1(ρ)(x))+J4, getting the following estimate for J4in a similar way |J4| ≤ rkρ0kL∞k∇ηkL∞(Br(x)) +r−1kρ0kL2. Thus, we conclude the following estimate for the factor L: |L(x)| ≤ ck∇ηkL∞(Br(x)) max j=1,2|Rj(ρ)|+ (rk∇ηkL∞(Br(x)) + 1)(k∇ηkL∞(Br(x))kρ0k∞+r−2kρ0k2). ANALYTICAL BEHAVIOR OF 2DPM 9 Using (3.10), we obtain by Gronwall’s lemma sup OT(Ω) |∇ρ(x, t)| ≤ sup Ω |∇ρ0|exp sup y∈ΩZT 0 M(t)dt, where M(t) is defined by M(t) = ck∇ηkL∞(Br(x)) max j=1,2|Rj(ρ)|+ (rk∇ηkL∞(Br(x)) + 1)(k∇ηkL∞(Br(x))kρ0k∞+r−2kρ0k2), with x= Φ(y, t). This concludes the proof of Theorem 3.4. Remark 3.6. The condition (3.15) depending on the Riesz transform is different that in 2DQG (see [8]). This appears because the integral kernels (1.9) in 2DPM are different to the Kernels in 2DQG. Now, we present a geometric conserved quantity that relates the curvature of the level sets and the magnitude |∇⊥ρ|in a similar way as in [6] (see references therein for more details). In particular, if we define the curvature of the level sets κby (3.18) κ(x, t) = (η· ∇η)·η⊥(x, t), where ηis the direction of ∇⊥ρ(see (3.12)), the following identity is satisfied (3.19) D(κ|∇⊥ρ|) Dt =∇⊥ρ· ∇β with (3.20) β(x, t) = (η· ∇v)·η⊥(x, t). Indeed, we now prove the identity (3.19). Since ∇⊥ρand |∇⊥ρ|satisfies (1.11) and (3.10) respectively, we get Dη Dt = (∇v)η− L η. Using (3.16), we obtain Dη Dt =β η⊥, with βdefined in (3.20). By the definition of κ(3.18) and the previous formula, we have Dκ Dt = ((∇v·η− L η)∇η)·η⊥+ (η·(β∇η⊥+η⊥⊗ ∇β− ∇η∇v)) ·η⊥−β(η· ∇η)·η and, simplifying, Dκ Dt = (∇β)η− L κ. Using this identity and (3.10), (3.19) is satisfied. Remark 3.7. The integral of the quantity κ|∇⊥ρ|over a region given by two different level sets is conserved along the time, i.e. (3.21) d dtZ {x:C1≤ρ(x,t)≤C2} κ|∇⊥ρ|dx= 0. This can be showed using the equation (3.19) and integrating by parts. Thus, in the case that |∇⊥ρ| is large by (3.21) the curvature κis small if the level sets do not oscillate. In all of our numerical experiments, we find no evidence of level set oscillations. On the contrary, we observe that the level sets are flattering where the gradient of ρis growing. 4. Numerical simulations Here we present two examples of numerical simulations for solutions of the 2-D mass balance in porous media with initial data in a period-cell [0,2π]2. Although periodic boundary conditions are rather unphysical, which does not matter because we are interested in the small-scale structures. 16 D. C ´ ORDOBA, F. GANCEDO AND R. ORIVE [11] D. C´ ordoba and C. Fefferman,Growth of solutions for QG and 2D Euler equations J. Amer. Math. Soc. 15 (2002), no. 3, 665–670. [12] W. E and C.-W. Shu,Small-scale structures in Boussinesq convection, Phys. Fluids 6 (1994), no. 1, 49–58. [13] T. Kato and Ponce,Commutator estimates and the Euler and Navier-Stokes equations, Communications on Pure and Applied Mathematics 41 (1998), 891–907. [14] H. Kozono and Y. Taniuchi,Limiting case of the Sobolev inequality in BMO, with application to the Euler equations, Communications in Mathematical Physics 214 (2000), 191–200. [15] A. J. Majda and A. L. Bertozzi,Vorticity and incompressible flow, Cambridge Univ. Press, Cambridge, 2002. [16] D.A. Nield and A. Bejan,Convection in porous media, Springer-Verlag, New York, 1999. [17] M. Price,Introducing Groundwater, 2nd ed. Chapman & Hall, London, 1996. [18] M. Siegel, R. E. Caflisch and S. Howison,Global existence, singular solutions, and ill-posedness for the Muskat problem, Comm. Pure Appl. Math. 57 (2004), no. 10, 1374–1411. [19] E.M. Stein,Harmonic Analysis, Princeton University Press, Princeton, NJ, 1993. [20] H. Vandeven,Family of spectral filters for discontinuous problems, J, Sci. Comput. 6 (1991), 159–192. Diego C´ ordoba and Francisco Gancedo Departamento de Matem´ aticas Instituto de Matem´ aticas y F´ ısica Fundamental Consejo Superior de Investigaciones Cient´ ıficas Serrano 123, 28006 Madrid, Spain. E-mail address:[email protected].es and [email protected] Rafael Orive Departamento de Matem´ aticas Facultad de Ciencias Universidad Aut´ onoma de Madrid Crta. Colmenar Viejo km. 15, 28049 Madrid, Spain. E-mail address:[email protected]