Analytical behavior of two-dimensional incompressible flow in porous media

JOURNAL OF MATHEMATICAL PHYSICS 48, 065206 共2007兲

Analytical behavior of two-dimensional incompressible flow in porous media Diego Córdobaa兲 and Francisco Gancedob兲 Departamento de Matemáticas, Instituto de Matemáticas y Física Fundamental, Consejo Superior de Investigaciones Científicas, Serrano 123, 28006 Madrid, Spain

Rafael Orivec兲 Departamento de Matemáticas, Facultad de Ciencias, Universidad Autónoma de Madrid, Crta. Colmenar Viejo km. 15, 28049 Madrid, Spain 共Received 19 July 2006; accepted 24 October 2006; published online 6 June 2007兲

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. © 2007 American Institute of Physics. 关DOI: 10.1063/1.2404593兴

I. 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 k is the matrix position-independent medium permeabilities in the different directions respectively divided by the viscosity, ␳ is the liquid density, g is 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 two-dimensional 共2D兲 fluid through porous media with two different constant densities and viscosities at each side of the interface is studied in Refs. 18 and 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

a兲

Electronic mail: [email protected] Electronic mail: [email protected] Electronic mail: [email protected]

b兲 c兲

0022-2488/2007/48共6兲/065206/19/$23.00

48, 065206-1

© 2007 American Institute of Physics

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-2

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

␾

冉

冊

D␳ ⳵␳ =␾ + v · ⵜ␳ = 0, Dt ⳵t

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␳ = 0, Dt

共1.1兲

v = − 共ⵜp + ␥␳兲.

共1.2兲

We close the system assuming incompressibility, i.e., 共1.3兲

div v = 0, therefore there exists a stream function ␺共x , t兲 such that

冉

冊

⳵␺ ⳵␺ , . ⳵x2 ⳵x1

v = ⵜ ⬜␺ ⬅ −

共1.4兲

Computing the curl of Eq. 共1.2兲, we get the Poisson equation for ␺, − ⌬␺ =

⳵␳ . ⳵x1

共1.5兲

A solution of this equation is given by the convolution of the Newtonian potential with the function ⳵x1␳,

␺共x,t兲 = −

1 2␲

冕

ln兩x − y兩

R2

⳵␳ 共y,t兲dy, ⳵y1

x 苸 R2 .

共1.6兲

Thus, the velocity v can be recovered from ␺ by the operator ⵜ⬜ by the two equivalent formulas, v共x,t兲 =

v共x,t兲 = PV

冕

R2

K共x − y兲ⵜ⬜␳共y,t兲dy,

x 苸 R2 ,

冕

1 H共x − y兲␳共y,t兲dy − 共0, ␳共x兲兲, 2 R2

共1.7兲

x 苸 R2 ,

共1.8兲

where the kernels K共·兲 and H共·兲 are defined by K共x兲 = −

1 x1 2␲ 兩x兩2

and

H共x兲 =

冉

冊

1 x 1x 2 x 2 − x 2 −2 4, 1 4 2 . 2␲ 兩x兩 兩x兩

共1.9兲

Differentiating Eq. 共1.1兲, we obtain the evolution equation for

冉

冊

共1.10兲

Dⵜ⬜␳ = 共ⵜv兲ⵜ⬜␳ . Dt

共1.11兲

ⵜ ⬜␳ ⬅ −

⳵␳ ⳵␳ , , ⳵x2 ⳵x1

which is given by

Taking the divergence of Eq. 共1.2兲, we get

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-3

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

− ⌬p =

⳵␳ ⳵x2

and the pressure can be obtained as in Eq. 共1.6兲. The objective of this work is to analyze the behavior of the solutions of the system 2DPM 关Eqs. 共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 Eq. 共1.10兲. A similar result is known in the three-dimensional Euler 共3D Euler兲 equation.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 equations7 and to the twodimensional quasigeostrophic equation 共2DQG兲.8 Applying these criteria, we find no evidence of formation of singularities in our numerical simulations. The paper is organized as follows. In Sec. II we study the analytical behavior of solutions with infinite energy. In Sec. III 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 Sec. IV we illustrate two numerical examples in which the analytical results are applied to show nonsingular solutions. II. SINGULARITIES WITH INFINITE ENERGY

Let the stream function ␺ be defined by

␺共x1,x2,t兲 = x2 f共x1,t兲 + g共x1,t兲.

共2.1兲

Note that under this hypothesis the solution of Eqs. 共1.1兲–共1.3兲 has infinite energy. We reduce Eqs. 共1.1兲–共1.3兲 to other system with respect to the functions f and g. From Eq. 共1.5兲 the density, apart from a constant, satisfies

␳共t,x1,x2兲 = − x2

⳵f ⳵g 共x1,t兲 − 共x1,t兲 = − x2 f x1 − gx1 , ⳵x1 ⳵x1

共2.2兲

and, by Eq. 共1.4兲, v verifies

冉

v共t,x1,x2兲 = − f共x1,t兲,x2

冊

⳵f ⳵g 共x1,t兲 + 共x1,t兲 = 共− f,x2 f x1 + gx1兲. ⳵x1 ⳵x1

共2.3兲

Therefore, the system 关Eqs. 共1.1兲–共1.3兲兴 under hypothesis 共2.1兲 is equivalent to 共f x兲t = f f xx − 共f x兲2 ,

共2.4兲

共gx兲t = fgxx − f xgx .

共2.5兲

共Here and in the sequel of the section, we denote with subscripts the derivatives with respect to x.兲 We note the nonlinear character of the first equation. Thus, our study of formation of singularities is concentrated in the solutions of Eq. 共2.4兲. The function g depends implicitly on f in Eq. 共2.5兲. Now, we show that the system 关Eqs. 共2.4兲 and 共2.5兲兴 is local well posed in the Sobolev spaces Hk0共0 , 1兲. Lemma 2.1: Let f 0 = f共x , 0兲 and g0 = g共x , 0兲 satisfy f 0x , g0x 苸 Hk0共0 , 1兲 with k 艌 1. Then, there exists T ⬎ 0 such that f x , gx 苸 C1共关0 , T兴 ; Hk0共0 , 1兲兲 are the unique solution of Eqs. (2.4) and (2.5). Proof: By Eq. 共2.4兲 and integrating by parts, we have

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-4

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

冕

1d 2 储f x储L2 = 2 dt

冕

1

f x f f xx −

0

3 2

冕

f xx f f xxx = −

1 2

1

f 3x = −

0

1 3

2

0

f 3x 艋 C储f x储L⬁储f x储L2 艋 C储f x储H1 . 0

Analogously, 1d 2 储f xx储L2 = − 2 dt

冕

冕

1 2 f xx fx −

0

1

0

冕

1 3

2

2 f xx f x 艋 C储f x储L⬁储f xx储L2 艋 C储f x储H1 . 0

0

We can repeat for all k 艌 1 and we obtain 1d 2 3 储f x储Hk 艋 C储f x储Hk . 0 0 2 dt Integrating in time, we get 储f 0x 储Hk

储f x储Hk 艋

0

1 − Ct储f 0x 储Hk

0

.

0

On the other hand, by Eq. 共2.5兲 and integrating by parts, we have for gx the following inequalities: 1d 2 储gx储L2 = 2 dt

冕

1

gxgxx f −

0

冕

1

3 2

g2x f x = −

0

冕

1 2

g2x f x 艋 储f x储L⬁储gx储H1 0

0

and 1d 2 储gxx储L2 = 2 dt

冕

冕

1

gxxgxxx f −

0

1

gxxgx f xx = −

0

1 2

冕

冕

1 2 gxx fx −

0

1 2

gxxgx f xx 艋 储f x储L⬁储gxx储L2

0

+ 储gxx储L2储gx储L⬁储f xx储L2 2

艋 储f x储H1储gx储H1 . 0

0

Thus, we obtain using Gronwall’s Lemma, 2

冉冕

2

储gx储Hk 艋 储g0x 储Hk exp C 0

0

冊

t

储f x储Hk ds , 0

0

and we have existence up to a time T = T共储f 0x 储Hk 兲. 0 In order to prove the uniqueness, let f x共x , t兲 = hx共x , t兲 − kx共x , t兲, with hx and kx two solutions of Eq. 共2.4兲 with the same initial data f 0x . Since hx and kx satisfy Eq. 共2.4兲 and integrating by parts, we have 1d 2 储f x储L2 = 2 dt =

冕 冕

1

f x共hhxx − kkxx兲 −

0 1

f xhf xx +

0

=−

1 2

冕

冕

1

0

1

f x fkxx −

0

1

0

冕

共f x兲2hx +

冕

f x共h2x − k2x 兲

冕

1

共f x兲2共hx + kx兲

0

1

0

f x fkxx −

冕

1

共f x兲2共hx + kx兲.

0

Thus, we get

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-5

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

1d 2 2 2 储f x储L2 艋 储kxx储L2储f x储L2储f储L⬁ + C共储hx储L⬁ + 储kx储L⬁兲储f x储L2 艋 C共储hx储H1 + 储kx储H1兲储f x储L2 0 0 2 dt and using Gronwall’s Lemma it follows hx = kx. Finally, we conclude the uniqueness of gx since Eq. 共2.5兲 is a linear differential equation. 䊏 The following result shows that the solution of Eq. 共2.4兲 blows up in finite time under certain conditions on the initial data. Proposition 2.2: Let f x be a solution of Eq. (2.4) with initial data satisfies f 0x 苸 H20共0 , 1兲 and minx f 0x ⬍ 0. Then, 储f x储L⬁ blows up in finite time T = −1 / minx f 0x . Proof: By the local existence result, we have f x 苸 C1共关0 , T兴 ; H2兲 傺 C1共关0 , T兴 ⫻ 关0 , 1兴兲. We consider the application m : 关0 , T兴 → R defined by m共t兲 = minx f x共x , t兲 = f x共xt , t兲. By the Rademacher Theorem, it follows that m is differentiable at almost every point. First, we calculate the derivative of m as in Refs. 5 and 9. Let s be a point of differentiability of m共t兲, then for ␶ ⬎ 0, m⬘共s兲 = lim

␶→0

m共s + ␶兲 − m共s兲 ␶

= lim

f x共xs+␶,s + ␶兲 − f x共xs,s兲 ␶

= lim

f x共xs+␶,s + ␶兲 − f x共xs,s + ␶兲 f x共xs,s + ␶兲 − f x共xs,s兲 + . ␶ ␶

␶→0

␶→0

Since f x共x , s + ␶兲 reaches the minimum at the point xs+␶, we obtain m⬘共s兲 艋 lim

␶→0

f x共xs,s + ␶兲 − f x共xs,s兲 = f xt共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 m⬘共s兲 = f xt共xs,s兲

almost everywhere.

We replace x for xs in Eq. 共2.4兲 and yields m⬘共s兲 = − f 2x 共xs,s兲 = − 共m共s兲兲2 , 䊏 due to f xx共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 f 0x 苸 H10 and assuming that

冕

1

f 0x 艋 0.

0

Thus, by Eq. (2.4), we have d dt

冕 冕 1

1

fx =

0

0

f f xx −

冕

1

共f x兲2 = − 2

0

冕

1

0

冉冕 冊 1

共f x兲2 艌 − 2

fx

2

.

0

Defining c共t兲 =

冕

fx,

and integrating, we get

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-6

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

c共t兲 艋

c共0兲 . 1 + 2tc共0兲

Then, c共t兲 blows up for c共0兲 ⬍ 0. In the case c共0兲 = 0, we have c⬘共t兲 ⬍ 0 for all t ⬎ 0, therefore, c共t兲 also blows up. Remark 2.4: Let x1 = xt be the point such that

䊏

f x共xt,t兲 = min f x共x,t兲, x

and consider x2 = 1 −

gx共xt,t兲 . f x共xt,t兲

Then, by Eq. (2.2), ␳共x1 , x2 , t兲 = −f x共xt , t兲 blows up in finite time by Proposition 2.2. Analogously, v defined in Eq. (2.3) blows up in finite time. 䊏 III. ANALYSIS OF 2DPM WITH FINITE ENERGY A. Local existence of 2DPM

We derive a reformulation of the system as an integrodifferential equation for the particle trajectories. Given a smooth field v共x , t兲, the particle trajectories ⌽共␣ , t兲 satisfy d⌽ 共␣,t兲 = v共⌽共␣,t兲,t兲, dt

⌽兩共␣,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 共Sec. 2.5 in Ref. 15兲 that Eq. 共1.11兲 implies the following formula: ⵜ⬜␳共⌽共␣,t兲,t兲 = ⵜ␣⌽共␣,t兲ⵜ⬜␳0共␣兲, where ⵜ⬜␳0 is 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 Eq. 共1.7兲 as v共⌽共␣,t兲,t兲 =

冕

R2

K共⌽共␣,t兲 − ⌽共␤,t兲兲ⵜ␣⌽共␤,t兲ⵜ⬜␳0共␤兲d␤ .

共3.2兲

We consider Eq. 共3.1兲 as an ordinary differential equation 共ODE兲 on a Banach space and using the 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 Sec. 4.1 in Ref. 15兲. In fact, we consider ⵜ⬜␳0 苸 C␦共R2兲, ␦ 苸 共0 , 1兲. Let B be the Banach space defined by B = 兵⌽:R2 → R2 such that 兩⌽共0兲兩 + 兩ⵜ␣⌽兩0 + 兩ⵜ␣⌽兩␦ ⬍ ⬁其, where 兩 · 兩0 is the L⬁ norm and 兩 · 兩␦ is the Hölder seminorm. Define OM , the open set of B, as

再

O M = ⌽ 苸 B兩 inf det ⵜ␣⌽共␣兲 ⬎ ␣苸R2

冎

1 and 兩⌽共0兲兩 + 兩ⵜ␣⌽兩0 + 兩ⵜ␣⌽兩␦ ⬍ M . 2

The mapping v共⌽兲, defined by Eq. 共3.2兲, satisfies the assumptions of the Picard Theorem, i.e., v is bounded and locally Lipschitz continuous on O M . As a consequence, for any M ⬎ 0 there exists T共M兲 ⬎ 0 and a unique solution

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-7

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

⌽ 苸 C1共共− T共M兲,T共M兲兲;O M 兲 to the particle trajectories 关Eqs. 共3.1兲 and 共3.2兲兴. Remark 3.1: The 2DPM [Eqs. (1.1)–(1.3)] has quantities conserved in time, the L p norm of ␳ for 1 艋 p 艋 ⬁, i.e., 储␳共t兲储 p = 储␳0储 p,

∀ t ⬎ 0,

1 艋 p 艋 ⬁.

共3.3兲

The velocity is obtained from ␳ by Eq. (1.8). These operators are singular integrals with CalderónZygmund kernels (see Ref. 19). Then for 1 ⬍ p ⬍ ⬁ the L p norm of the velocity is bounded for any time t ⬎ 0. B. 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 Chap. IV in Ref. 19 for an introduction of this function space兲. Theorem 3.2: Let ␳ be the solution of Eqs. (1.1)–(1.3) with initial data ␳0 苸 Hs共R2兲 with s ⬎ 2. Then, the following are equivalent: 共a兲 共b兲

The interval 关0 , ⬁兲 is the maximal interval of Hs existence for ␳. The quantity

冕

T

储ⵜ␳储BMO共t兲dt ⬍ ⬁,

∀ T ⬎ 0.

共3.4兲

0

Proof: We denote the operator ⌳s by ⌳s ⬅ 共−⌬兲s/2. Since the fluid is incompressible, we have for s ⬎ 2 1d s 2 储⌳ ␳储L2 = − 2 dt

冕

R2

⌳s␳⌳s共v ⵜ ␳兲dx = −

冕

R2

⌳s␳共⌳s共v ⵜ ␳兲 − v⌳s共ⵜ␳兲兲dx,

艋 C储⌳s␳储L2储⌳s共v ⵜ ␳兲 − v⌳s共ⵜ␳兲储L2 . Using the estimate 共see Ref. 13兲, 储⌳s共fg兲 − f⌳s共g兲储Lp 艋 c共储ⵜf储L⬁储⌳s−1g储Lp + 储⌳s f储Lp储g储L⬁兲,

1 ⬍ p ⬍ ⬁,

we obtain for p = 2, 1d s 2 2 储⌳ ␳储L2 艋 C共储ⵜv储L⬁ + 储ⵜ␳储L⬁兲储⌳s␳储L2 . 2 dt Integrating, we get for any t 艋 T,

冉冕

T

储⌳s␳储L2 艋 储⌳s␳0储L2 exp C

共3.5兲

冊

共储ⵜv储L⬁ + 储ⵜ␳储L⬁兲 .

0

共3.6兲

Now, we use the following inequality given in Ref. 14. Let f 苸 Ws,p with 1 ⬍ p ⬍ ⬁ and s ⬎ 2 / p, then, there exists a constant C = C共p , s兲 such that 储f储L⬁ 艋 C共1 + 储f储BMO共1 + ln+储f储Ws,p兲兲,

共3.7兲

where ln+共x兲 = max共0 , ln共x兲兲. Therefore, for s ⬎ 2 we have 储ⵜ␳储L⬁ 艋 C共1 + 储ⵜ␳储BMO共1 + ln+储ⵜ␳储Hs−1兲兲, and from Eq. 共3.6兲, we obtain that

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-8

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

储ⵜ␳储L⬁ 艋 C共1 + ln+共储␳0储Hs兲储ⵜ␳储BMO

冕

T

共储ⵜv储L⬁ + 储ⵜ␳储L⬁兲dt兲.

共3.8兲

0

On the other hand, applying Eq. 共3.7兲 for v 苸 Hs共R2兲, we have 储ⵜv储L⬁ 艋 C共1 + 储ⵜv储BMO共1 + ln+储ⵜv储Hs−1兲兲. Since v satisfies Eq. 共1.8兲 and the singular integrals are bounded operators in BMO 共see Ref. 19兲, we get 储ⵜv储L⬁ 艋 C共1 + 储ⵜ␳储BMO共1 + ln+储ⵜ␳储Hs−1兲兲, and, using Eq. 共3.6兲, we obtain

冉

储ⵜv储L⬁ 艋 C 1 + ln+共储␳0储Hs兲储ⵜ␳储BMO

冕

冊

T

共储ⵜv储L⬁ + 储ⵜ␳储L⬁兲dt .

0

共3.9兲

From Eqs. 共3.8兲 and 共3.9兲 follows

冉

储ⵜv储L⬁ + 储ⵜ␳储L⬁ 艋 C 1 + ln+共储␳0储Hs兲储ⵜ␳储BMO

冕

冊

T

共储ⵜv储L⬁ + 储ⵜ␳储L⬁兲dt .

0

Applying Gronwall’s inequality, we have

冕

冉

T

共储ⵜv储L⬁ + 储ⵜ␳储L⬁兲dt 艋 CT exp ln+共储␳0储Hs兲

0

冕

T

0

冊

储ⵜ␳储BMOdt ,

and so 共a兲 is a consequence of 共b兲. Finally, due to the inequality 储ⵜ␳储BMO 艋 储ⵜ␳储H1 , we conclude that 共a兲 implies 共b兲. Remark 3.3: Using that

䊏 储ⵜ␳储BMO 艋 C储ⵜ␳储L⬁ ,

we get a blow-up characterization for numerical simulations. C. Geometric constraints on singular solutions

From Eq. 共1.1兲 it follows that the level sets, ␳ = constant, move with the fluid flow. Then ⵜ⬜␳, defined in Eq. 共1.10兲, is tangent to these level sets. For the 2DPM, the infinitesimal length of a level set for ␳ is given by 兩ⵜ⬜␳兩 and from Eq. 共1.11兲, the evolution equation for the infinitesimal arc length is given by D兩ⵜ⬜␳兩 = L兩ⵜ⬜␳兩. Dt The factor L共x , t兲 is defined through by L共x,t兲 =

再

D␩ · ␩, 0,

␩⫽0

␩ = 0,

共3.10兲

冎

共3.11兲

where the direction of ⵜ⬜␳ is denoted by

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-9

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

␩=

ⵜ ⬜␳ , 兩ⵜ⬜␳兩

共3.12兲

and D共x , t兲 is the symmetric part of the deformation matrix defined by D = 共Dij兲 =

冋冉

1 ⳵vi ⳵v j + 2 ⳵x j ⳵xi

冊册

共3.13兲

.

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 Euler7 and to 2DQG.8 Recall that ␩ is the direction field tangent to the level sets of ␳ defined 关Eq. 共3.12兲兴. Analogously to Ref. 8 a set ⍀ is smoothly directed if there exists ␦ ⬎ 0 such that sup ¯ x苸⍀

冕

T

0

储ⵜ␩共·,t兲储L⬁共B 2

␦共⌽共x,t兲兲兲

dt ⬍ ⬁,

共3.14兲

where B␦共x兲 = 兵y 苸 R2: 兩x − y兩 ⬍ ␦其,

¯ = 兵x 苸 ⍀;兩ⵜ␳ 共x兲兩 ⫽ 0其, ⍀ 0

and ⌽ is the particle trajectories map. We define ⍀共t兲 = ⌽共⍀ , t兲 and OT共⍀兲 the semiorbit, i.e., OT共⍀兲 = 艛 兵t其 ⫻ ⍀共t兲. 0艋t艋T

Theorem 3.4: If ⍀ is smoothly directed and

冕

T

储R j␳储L⬁共t兲dt ⬍ ⬁,

j = 1,2,

∀ T ⬎ 0,

共3.15兲

0

where R j denotes the Riesz transform in the direction x j, then sup 兩ⵜ␳共x,t兲兩 ⬍ ⬁.

OT共⍀兲

Remark 3.5: Using 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 L defined in Eq. 共3.11兲. We start by computing the full gradient of the velocity v. From formula 共1.7兲, v共x兲 =

冕

R2

K共y兲ⵜ⬜␳共x − y兲dy,

we have ⵜv共x兲 =

冕

R2

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 ⵜv共x兲 =

1 PV 2␲

冕

R2

共ⵜy␳共x − y兲兲 丢 ˜y

冉

冊

0 0 dy 1 − , 兩y兩2 2 共⳵␳/⳵x 兲共x兲 共⳵␳/⳵x 兲共x兲 1 2

共3.16兲

where ˜y is the unit vector defined by

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-10

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

冉

˜y = −

冊

2y 1y 2 y 21 − y 22 , . 兩y兩2 兩y兩2

By definition of ␩ in Eq. 共3.12兲, we have ␩ · ⵜ␳ = 0. Thus, computing we get 1 PV 2␲

L共x兲 =

冕

R2

˜ · ␩共x兲兲兵共␩共x − y兲 · ␩⬜共x兲兲其兩ⵜ⬜␳共x − y兲兩 共y

dy . 兩y兩2

共3.17兲

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␲

冏冕

R2

˜ · ␩共x兲兲共ⵜ⬜␳共x − y兲 · ␩⬜共x兲兲 共1 − ␾共兩y兩2/r2兲兲共y

冏

dy . 兩y兩2

Integrating by parts and using the Cauchy-Schwartz inequality, we get I1 艋

C C 2 储 ␳ 储 L2 艋 2 储 ␳ 0储 L2 . r r

We have for any 兩y 兩 ⬍ r, 兩␩共x − y兲 · ␩⬜共x兲兩 艋 兩y兩储ⵜ␩储L⬁共Br共x兲兲 . Applying this in the integral I2, we get 兩I2兩 艋

冕

兩ⵜ⬜␳共x − y兲兩␾共兩y兩2/r2兲

dy 储ⵜ␩储L⬁共Br共x兲兲 . 兩y兩

We integrate by parts and decompose,

冕

兩ⵜ⬜␳共x − y兲兩␾共兩y兩2/r2兲

dy = 兩y兩

冕

冉

␳共x − y兲ⵜ⬜ ␩共x − y兲␾共兩y兩2/r2兲

冊

1 dy = J1 + J2 + J3 , 兩y兩

where J1 =

J2 =

J3 =

冕

␳共x − y兲␩共x − y兲ⵜ⬜共␾共兩y兩2/r2兲兲

dy , 兩y兩

冕

␳共x − y兲ⵜ⬜共␩共x − y兲兲␾共兩y兩2/r2兲

dy , 兩y兩

冕

␳共x − y兲␩共x − y兲␾共兩y兩2/r2兲

共− y 2,y 1兲 dy, 兩y兩3

obtaining the following estimates: 兩J1兩 艋 c储␳0储L⬁

and

兩J2兩 艋 cr储␳0储L⬁储ⵜ␩储L⬁共Br共x兲兲 .

The J3 term can be bounded using the identity

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-11

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

J3 = ␩共x兲共− R2共␳兲共x兲,R1共␳兲共x兲兲 + J4 , getting the following estimate for J4 in a similar way: 兩J4兩 艋 r储␳0储L⬁储ⵜ␩储L⬁共Br共x兲兲 + r−1储␳0储L2 . Thus, we conclude the following estimate for the factor L: 兩L共x兲兩 艋 c关储ⵜ␩储L⬁共Br共x兲兲max兩R j共␳兲兩 + 共r储ⵜ␩储L⬁共Br共x兲兲 + 1兲共储ⵜ␩储L⬁共Br共x兲兲储␳0储⬁ + r−2储␳0储2兲兴. j=1,2

Using Eq. 共3.10兲, we obtain by Gronwall’s Lemma,

冉 冕

T

sup 兩ⵜ␳共x,t兲兩 艋 sup兩ⵜ␳0兩exp sup

OT共⍀兲

⍀

y苸⍀

0

冊

M共t兲dt ,

where M共t兲 is defined by M共t兲 = c关储ⵜ␩储L⬁共Br共x兲兲max兩R j共␳兲兩 + 共r储ⵜ␩储L⬁共Br共x兲兲 + 1兲共储ⵜ␩储L⬁共Br共x兲兲储␳0储⬁ + r−2储␳0储2兲兴, j=1,2

with x = ⌽共y , t兲. This concludes the proof of Theorem 3.4. 䊏 Remark 3.6: Condition (3.15) depending on the Riesz transform is different from that in 2DQG (see Ref. 8). This appears because the integral kernels [Eq. (1.9)] in 2DPM are different from 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 Ref. 6 共see references therein for more details兲. In particular, if we define the curvature of the level sets ␬ by

␬共x,t兲 = 共␩ · ⵜ␩兲 · ␩⬜共x,t兲,

共3.18兲

where ␩ is the direction of ⵜ⬜␳ 关see Eq. 共3.12兲兴, the following identity is satisfied: D共␬兩ⵜ⬜␳兩兲 = ⵜ ⬜␳ · ⵜ ␤ , Dt

共3.19兲

␤共x,t兲 = 共␩ · ⵜv兲 · ␩⬜共x,t兲.

共3.20兲

with

Indeed, we now prove the identity 共3.19兲. Since ⵜ⬜␳ and 兩ⵜ⬜␳兩 satisfies Eqs. 共1.11兲 and 共3.10兲, respectively, we get D␩ = 共ⵜv兲␩ − L␩ . Dt Using Eq. 共3.16兲, we obtain D␩ = ␤␩⬜ , Dt with ␤ defined in Eq. 共3.20兲. By the definition of ␬ 关Eq. 共3.18兲兴 and the previous formula, we have D␬ = 共共ⵜv · ␩ − L␩兲 ⵜ ␩兲 · ␩⬜ + 共␩ · 共␤ ⵜ ␩⬜ + ␩⬜ 丢 ⵜ␤ − ⵜ␩ ⵜ v兲兲 · ␩⬜ − ␤共␩ · ⵜ␩兲 · ␩ Dt and simplifying,

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-12

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

FIG. 1. Evolution of the density in case 1 for times t = 0, 3, 6, and 8.5.

D␬ = 共ⵜ␤兲␩ − L␬ . Dt Using this identity and Eq. 共3.10兲, Eq. 共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., d dt

冉冕

兵x:C1艋␳共x,t兲艋C2其

冊

␬兩ⵜ⬜␳兩dx = 0.

共3.21兲

This can be shown using Eq. (3.19) and integrating by parts. Thus, in the case that 兩ⵜ⬜␳兩 is large by Eq. (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. IV. NUMERICAL SIMULATIONS

Here we present two examples of numerical simulations for solutions of the 2D 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. The numerical method. We solve Eqs. 共1.1兲–共1.3兲 numerically on a 2␲-periodic cell with a spectral method with smoothing. This numerical method is similar to the scheme developed by E and Shu12 for incompressible flows and also used for the quasigeostrophic active scalar in Ref. 8. This algorithm is the standard Fourier-collocation method 共see Ref. 4兲 with smoothing. Roughly speaking, the differentiation operator is approximated in the Fourier space, while the nonlinear operations such as v · ⵜ␳ are done in the physical space.

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-13

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

FIG. 2. Evolution of around contour ␳共x1 , x2兲 = 1 in case 2 for times t = 0, 1.5, 3, 4.5, 6, and 7.5.

We smooth the gradients adding filters to the spectral method in order for the numerical solutions do not degrade catastrophically. A way of adding the filters is to replace the Fourier multiplier ik j by ik j␸共兩k j兩兲, where

␸共k兲 = e−a共k/N兲

b

for 兩k兩 艋 N.

Here N is the numerical cutoff for the Fourier modes and a and b are chosen with the machine accuracy 共see Ref. 20兲. For the temporal discretization, we use Runge-Kutta methods of various orders. In our case, we have no explicit temporal dependence and we get a Runge-Kutta method of order of 4 that requires at most three levels of storage 共see in Ref. 4, p. 109兲. We present the numerical approximation with an initial resolution of 共256兲2 Fourier modes. We refine this resolution when the growth of 储ⵜ␳储L⬁ is substantial to give additional insight preserving the relation space-time. We conclude our numerical simulations with a resolution of 共8192兲2 Fourier modes.

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-14

Córdoba, Gancedo, and Orive

J. Math. Phys. 48, 065206 共2007兲

FIG. 3. Evolution of the level sets −0.999, −0.99 共on the left of the figure兲 and 0.99, 0.999 共on the right兲 of the density in case 1 for times t = 0, 3, and 6.

The computational part of this work was performed on HPC320 共cluster of eight servers with 32 processors Alpha EV68 1 GHz兲 of the Centro de Supercomputación de Galicia 共CESGA兲. We used the MATLAB routines to obtain the calculations. Case 1. We consider the initial datum,

␳共x1,x2兲 = sin共x1兲sin共x2兲. The time step is ⌬t = 0.025 from t = 0 to t = 4.0 with a = 4.5 and b = 2.3, stopping the experiment with ⌬t = 0.001 25 and a = 9.1 and b = 7.1. During the simulation the ratio h / ⌬t is preserved getting a finer resolution as the gradients are growing. In this way the method conserves in time the L⬁ norm of the density. Figure 1 presents the density at times t = 0, 3, 6, and 8.5 with a numerical resolution of 共256兲2, 共256兲2, 共1024兲2, and 共4096兲2, respectively. The initial data have a hyperbolic saddle point that does not present a nonlinear behavior as time evolves. In all our numerical simulations we do not find

FIG. 4. Two level sets that are approaching each other.

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-15

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

FIG. 5. Evolution of the L⬁ norms of ⵜ␳, the velocity v, and the Riesz transforms 共R1␳ , R2␳兲 for case 1.

any saddle point structures that present stronger front formation than in case 2 共see Fig. 2兲. Nevertheless, case 1 shows a front formation in the elliptic plots and approaching the maximum with the minimum 共see Fig. 3兲 as time evolves. If we choose two level sets that are approaching each other we have the following scenario where ␦ = ␦共x2 , t兲 is the distance between the two counters 共see Fig. 4兲. From previous work 共see Ref. 10兲 it is easy to check that in order for the two graphs f l and f r to collapse at time T in any interval x2 苸 关a , b兴, i.e., lim 关f r共x2,t兲 − f l共x2,t兲兴 = 0,

t→T−

∀ x2 苸 关a,b兴,

it is necessary that

冕

T

储v储⬁共s兲ds = ⬁.

共4.1兲

0

In Fig. 5 we see no evidence for quantity 兰T0 储v储⬁共s兲ds to blow up in finite time. Nevertheless, we can obtain an estimate of how close the two graphs approach each other without any assump-

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-16

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

FIG. 6. Evolution of the density in case 2 for times t = 0, 3, 6, and 9.

tion on the velocity. From Fig. 3 it seems reasonable to assume that the minimum and maximum of ␦ are comparable that means there exists and constant c ⬎ 0 such that min ␦共x2,t兲 艋 c max ␦共x2,t兲,

∀ x2 苸 关a,b兴.

Then we can obtain an evolution equation for the area A = A共t兲,

A共t兲 =

1 b−a

冕

b

关f r共x2,t兲 − f l共x2,t兲兴dx2 ,

a

in between the two graphs that satisfies 共see Ref. 11 for more details兲

冏 冏

dA C 共t兲 艋 sup 兩␺共f r共x2,t兲,x2,t兲 − ␺共f l共x2,t兲,x2,t兲兩, dt b − a a艋x2艋b

共4.2兲

where ␺ is the stream function 关Eq. 共1.6兲兴,

␺共x,t兲 = −

1 2␲

冕

R2

x1 − z1 ␳共z,t兲dz, 兩x − z兩2

x 苸 R2 .

Using this formula we obtain that for any x , y 苸 R2, 兩␺共x,t兲 − ␺共y,t兲兩 艋 C共储␳0储⬁,储␳0储L2兲兩x − y兩共1 − ln兩x − y兩兲,

共4.3兲

due to the following estimates:

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-17

J. Math. Phys. 48, 065206 共2007兲

2D incompressible flow in porous media

FIG. 7. Evolution of L⬁ norms of ⵜ␳, the velocity v, and the Riesz transforms 共R1␳ , R2␳兲 for case 2.

兩␺共x,t兲 − ␺共y,t兲兩 =

冊 冏 冏 冕冉 冏 冕 冏冏 冕 冏冏 冕 冏

艋

1 2␲

1 2␲

R2

x1 − z1 y 1 − z1 − ␳共z,t兲dz 兩x − z兩2 兩y − z兩2 +

B2r共x兲

1 2␲

+

B2共x兲−B2r共x兲

1 2␲

B2共x兲

= I1 + I2 + I3 ,

where r = 兩x − y兩 and I1 艋 C储␳0储⬁兩x − y兩, I2 艋 C储␳0储⬁兩x − y兩

冕

2

s−1ds 艋 C储␳0储⬁兩x − y兩共− ln兩x − y兩兲,

2r

I3 艋 C储␳0储L2兩x − y兩. Then, using Eq. 共4.3兲 in Eq. 共4.2兲 we get that the area A共t兲 is bounded by

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-18

J. Math. Phys. 48, 065206 共2007兲

Córdoba, Gancedo, and Orive

FIG. 8. Level contour 共0.98,0.99,1.0,1.01,1.02兲 of ␳ in case 2 for times t = 7.5, 8, 8.5, and 9.

t

A共t兲 艌 A0e−Ce . In order to apply the global solution criterion 共Theorem 3.2兲, in Fig. 5 is plotted the logarithm of the L⬁ norms of ⳵1␳ and ⳵2␳ showing an exponential growth. The hypothesis 共3.15兲 of Theorem 3.4 is checked in this example 共see Fig. 5兲, so that using this result no singularity is possible due to the variation of the direction field tangent to the level sets is smooth. Case 2. In this case, the initial datum is

␳共x1,x2兲 = sin共x1兲cos共x2兲 + cos共x1兲. The time step is ⌬t = 0.025 from t = 0 to t = 4.5 with a = 4.5 and b = 2.3, stopping the experiment with ⌬t = 0.001 and a = 11.4 and b = 8. The L⬁ norm is preserved during the simulation of this case. Figure 6 presents the density at times t = 0, 3, 6, and 9 with a numerical resolution of 共256兲2, 共256兲2, 共1024兲2, and 共8192兲2, respectively. In Fig. 7, we show log plots of max兩⳵x1␳兩 and max兩⳵x2␳兩, where we obtain an analogous exponential growth as in case 1. This growth is not sufficient to guarantee a singularity. The initial data for the density scalar clearly have a hyperbolic saddle and the numerical solution develops a front as time evolves 共see Figs. 2 and 8兲. We observe that the front does not develop nonlinear or potentially singular structure as time evolves. Where ␩ is smoothly directed we observe the highest growth of 储ⵜ␳储L⬁. Where ␩ changes rapidly we obtain less growth of 储ⵜ␳储L⬁. Using Theorem 3.4 we show no evidence of singularities. ACKNOWLEDGMENTS

The authors were partially supported by Grant Nos. MTM2005-05980 of the MEC 共Spain兲 and S-0505/ESP/0158 of the CAM 共Spain兲. Two of the authors 共D.C. and F.G.兲 were partially

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

065206-19

2D incompressible flow in porous media

J. Math. Phys. 48, 065206 共2007兲

supported by Grant No. PAC-05-005-2 of the JCLM 共Spain兲. Another author 共R.O.兲 was partially supported by Grant No. MTM2005-00714 of the MEC 共Spain兲. 1

D. Ambrose, “Well-posedness of two-phase Hele-Shaw flow without surface tension,” Eur. J. Appl. Math. 15, 597–607 共2004兲. J. T. Beale, T. Kato, and A. Majda, “Remarks on the breakdown of smooth solutions for the 3-D Euler equations,” Commun. Math. Phys. 94, 61–66 共1984兲. 3 J. Bear, Dynamics of Fluids in Porous Media 共Elsevier, New York, 1972兲. 4 C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics 共Springer-Verlag, New York, 1988兲. 5 A. Constantin, and J. Escher, “Wave breaking for nonlinear nonlocal shallow water equations,” Acta Math. 181, 229–243 共1998兲. 6 P. Constantin, “Geometric statistics in turbulence,” SIAM Rev. 36, 73–98 共1994兲. 7 P. Constantin, C. Fefferman, and A. Majda, “Geometric constraints on potential singularity formulation in the 3-D Euler equations,” Commun. Partial Differ. Equ. 21, 559–571 共1996兲. 8 P. Constantin, A. Majda, and E. Tabak, “Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar,” Nonlinearity 7, 1495–1533 共1994兲. 9 A. Córdoba, and D. Córdoba, “A maximum principle applied to quasi-geostrophic equations,” Commun. Math. Phys. 249, 511–528 共2004兲. 10 D. Córdoba, and C. Fefferman, “Scalars convected by a two-dimensional incompressible flow,” Commun. Pure Appl. Math. 55, 255–260 共2002兲. 11 D. Córdoba, and C. Fefferman, “Growth of solutions for QG and 2D Euler equations,” J. Am. Math. Soc. 15, 665–670 共2002兲. 12 W. E, and C.-W. Shu, “Small-scale structures in Boussinesq convection,” Phys. Fluids 6, 49–58 共1994兲. 13 T. Kato, and G. Ponce, “Commutator estimates and the Euler and Navier-Stokes equations,” Commun. Pure Appl. Math. 41, 891–907 共1998兲. 14 H. Kozono, and Y. Taniuchi, “Limiting case of the Sobolev inequality in BMO, with application to the Euler equations,” Commun. Math. Phys. 214, 191–200 共2000兲. 15 A. J. Majda, and A. L. Bertozzi, Vorticity and Incompressible Flow 共Cambridge University 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 and Hall, London, 1996兲. 18 M. Siegel, R. E. Caflisch, and S. Howison, “Global existence, singular solutions, and ill-posedness for the Muskat problem,” Commun. Pure Appl. Math. 57, 1374–1411 共2004兲. 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, 159–192 共1991兲. 2

Downloaded 10 Nov 2008 to 161.111.22.69. Redistribution subject to ASCE license or copyright; see http://jmp.aip.org/jmp/copyright.jsp