Application of variational inequalities to the moving-boundary problem in a fluid model for biofilm growth

Nonlinear Analysis 70 (2009) 3658–3664

Contents lists available at ScienceDirect

Nonlinear Analysis journal homepage: www.elsevier.com/locate/na

Application of variational inequalities to the moving-boundary problem in a fluid model for biofilm growthI Niels Chr. Overgaard ∗ Applied Mathematics Group, School of Technology, Malmö University, SE-205 06 Malmö, Sweden

article

info

Article history: Received 22 May 2006 Accepted 29 July 2008 MSC: 35R35 49J10 92B05 Keywords: Biofilm model Moving-boundary problem Variational inequality

a b s t r a c t We consider a moving-boundary problem associated with the fluid model for biofilm growth proposed by J. Dockery and I. Klapper, Finger formation in biofilm layers, SIAM J. Appl. Math. 62 (3) (2001) 853–869. Notions of classical, weak, and variational solutions for this problem are introduced. Classical solutions with radial symmetry are constructed, and estimates for their growth given. Using a weighted Baiocchi transform, the problem is reformulated as a family of variational inequalities, allowing us to show that, for any initial biofilm configuration at time t = 0 (any bounded open set), there exists a unique weak solution defined for all t ≥ 0. © 2008 Elsevier Ltd. All rights reserved.

1. Introduction Given a bounded open set B ⊂ Rn , n > 0 fixed, let pB denote the solution of the boundary value problem

(−∆ + 1)pB = 1 in B, and pB = 0 on ∂ B.

(1)

The purpose of this paper is to show that the following moving-boundary problem admits solutions in a certain weak sense. (P) Given an bounded open set B ⊂ Rn and a number T > 0, find a set-valued function [0, T ] 3 t 7→ Bt ⊂ Rn such that (i) Bt is bounded and open for any 0 ≤ t ≤ T , (ii) B0 = B, and (iii) if pt = pBt is the solution of (1), then the free boundary Γt = ∂ Bt moves with the normal velocity

Γ˙ t = −

∂ pt , ∂n

(2)

where ∂/∂ n denotes the outward normal derivative on Γt . This problem has its origin in mathematical (micro-)biology, namely the biofilm model proposed by Dockery and Klapper [4]. The model describes the expansion of a bacterial colony caused by cell divisions. Mathematically, the problem may be viewed as a generalization of the Hele-Shaw squeeze film model (cf. [11] or [8, p. 2561]) in which the Laplacian is replaced by the modified Helmholtz operator. Our analysis is based on a reformulation of (P) as a (time-indexed) sequence of variational inequalities. This idea was pioneered by Elliott and Janovský [5] and further developed by Gustafsson [6] in the case of moving-boundary problems for Hele-Shaw flows. Our exposition is based on suitable modifications of ideas found in Gustafsson’s paper.

I Work supported by the Swedish Knowledge Foundation: http://www.kks.se/.

∗

Tel.: +46 40 665 71 32; fax: +46 665 76 46. E-mail address: [email protected].

0362-546X/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.na.2008.07.021

N.Chr. Overgaard / Nonlinear Analysis 70 (2009) 3658–3664

3659

2. The fluid model of Dockery and Klapper We recall the Dockery–Klapper model and explain how Problem (P) is derived. Let Bt denote the domain occupied by the biofilm at time t. A nutrient, e.g. oxygen, is diffusing into the biofilm from the surrounding bulk fluid (Bt ’s complement), where it is consumed at a rate proportional to its concentration S. In suitable dimension-free units, the governing equation becomes,

∂S − ∆S = −S in Bt , and S = 1 on Γt := ∂ Bt , ∂t

(3)

where the second equation expresses the assumption that the nutrient concentration is constant, S = 1, in the bulk fluid. As the bacteria multiply by cell division, a fixed fraction (here = 1) of the consumed nutrient is transformed into new biomass. If the biomass is regarded as a viscous fluid of fixed density, then the only way to accommodate the new cells is to let the biofilm expand. The expansion is described by a volumetric flow field u generated by a pressure p according to Darcy’s law, u = −∇ p. We set p = 0 in the bulk fluid. Now, conservation of mass states that div u = added biomass per unit volume = S , hence the pressure must satisfy the boundary value problem:

− ∆p − S = 0 in Bt ,

and p = 0

on Γt .

(4)

Based on the observation that the diffusion process is much faster than the growth process, cf. e.g. [10, p. 401], the nutrient field is taken to be the steady state solution of (3). By adding the pressure equation (4) to the nutrient equation (3), with ∂ S /∂ t = 0, we find that

−∆(p + S ) = 0 in Bt and p + S = 1 on Γt , hence p + S = 1 throughout Bt . Thus, S in (4) can be eliminated, resulting in Eq. (1) for the pressure. Finally, the normal velocity of the biofilm interface is obtained as the normal component of the flow, Γ˙ t = n · u = −∂ p/∂ n on Γt , where n denotes the outward unit normal of Bt . This concludes the derivation of the moving-boundary problem (P). 3. Classical solutions with radial solution In this section, a classical formulation of the moving-boundary problem of (P) is given, and radially symmetric classical solutions constructed. For a given number T > 0, let {Bt }0≤t ≤T denote a set-valued mapping [0, T ] 3 t 7→ Bt ⊂ Rn , where Bt is a bounded open set for each t ∈ [0, T ]. If Bt has a C 2 -boundary, and the pressure pBt is defined by (1), then the maximum principle implies ∂ pBt /∂ n < 0 on ∂ B, hence Bt will be expanding by the moving-boundary condition (2). This motivates the following definition: {Bt }0≤t ≤T is a classical solution to (P) if there exists a C 2 -function f = f (x, t ) : U × [0, T ], for some open set U with BT ⊂ U ⊂ Rn , which satisfies: (i) Bt = {x ∈ Rn : f (x, t ) > 0} for each t ∈ [0, T ]. (ii) Zero is a regular value of f (·, t ) : U → R for each t ∈ [0, T ], i.e., ∇x f 6= 0 on Γt . In particular Γt = {x ∈ Rn : f (x, t ) = 0} is of class C 2 . (iii) If p(x, t ) = pBt (x) is the solution of (1), then

∂f = ∇x p · ∇x f ∂t

on Γt for all t ∈ [0, T ].

(5)

Since the normal velocity of the moving boundary t 7→ Γt is related to f by

Γ˙ t =

∂ f /∂ t |∇x f |

on Γt ,

(6)

and n = ∇x f /|∇x f |, (5) is, of course, nothing but the moving-boundary condition (2). To construct radially symmetric classical solutions, we assume the initial set is B = {x ∈ Rn : |x| < R0 } and seek a non-negative function R : [0, ∞) → R with R(0) = R0 , such that Bt = {x : |x| < R(t )}, 0 ≤ t ≤ T , solves (P). First, we determine the pressure (1) for any ball B = {x : |x| < R}. If we set pB (x) = gR (r ), r = |x|, then gR has to solve the following boundary value problem on 0 ≤ r ≤ R: 1 r n−1

[r n−1 gR0 (r )]0 + gR (r ) = 1,

gR0 (0) = 0 and gR (R) = 0.

(7)

The solution can be expressed in terms of Bessel functions of the first kind, Jν (z ). In fact, if we set Gν (z ) = (iz /2)−ν Jν (iz ),

(8)

3660

N.Chr. Overgaard / Nonlinear Analysis 70 (2009) 3658–3664

then it is easy to check that the solution of (7) is gR (r ) = 1 −

Gν (r )

(0 ≤ r ≤ R; ν = (n − 2)/2).

Gν (R)

(9)

Since Γ˙ t = R0 (t ) and ∂ p/∂ n = gR0 (R), we get a classical solution of (P) if R0 (t ) =

G0ν (R(t ))

(t > 0; ν = (n − 2)/2), and R(0) = R0 . (10) √ For n = 1, G−1/2 (r ) = (cosh r )/ π , and the initial value problem (10) has the explicit solution R(t ) = t + ln[sinh R0 + (e−2t + sinh2 R0 )1/2 ]. Note that R(t ) ≤ t + R0 , for all t ≥ 0, and R0 (t ) is monotonically increasing with R0 (t ) → 1 as t → ∞. Gν (R(t ))

,

The latter means that the biofilm interface approaches a limiting speed monotonically. These observations are true for radial solutions in all dimensions n ≥ 1. This is a consequence of the following result, whose proof is given in the Appendix: Proposition 1. For ν ≥ − 12 the function G0ν (r )/Gν (r ), r ∈ R, is strictly monotonically increasing, G0ν (0)/Gν (0) = 0, and limr →±∞ G0ν (r )/Gν (r ) = ±1. The following theorem, which is used in Section 5, summarizes these results: Theorem 2. For t ≥ 0 set Bt = {x : |x| < R(t )}, where R(t ) is the solution of (10). Then, for any T > 0, the set-valued map {Bt }0≤t ≤T is a classical solution of (P) with initial data B0 = {x : |x| < R0 }. (Take e.g. f (x, t ) = R(t )2 − |x|2 .) This solution satisfies the bound Bt ⊂ {x : |x| < t + R0 }, for all t ≥ 0, and the normal velocity Γ˙ t of the moving boundary increases monotonically to the limiting speed 1 as the diameter of the biofilm tends to infinity. 4. Weak solutions Following Gustafsson [6] we now introduce a type of weak solution of (P), which is more flexible than the classical solutions of the previous section. The idea is to replace the sets in a classical solution, {Bt }0≤t ≤T , by characteristic functions:

χt (x) = χBt (x) =



1 0

if x ∈ Bt otherwise.

We can always find a bounded open set U with smooth boundary, for instance a large open ball, such that the closure of Bt is contained in U for all t ∈ [0, T ]. Clearly χt ∈ L2 (U ), and in the following {χt }0≤t ≤T is used as a short-hand for the mapping [0, T ] 3 t 7→ χt ∈ L2 (U ). Definition 3. A mapping {χt }0≤t ≤T is called a weak solution of (P) if the function [0, T ] 3 t 7→ ut ∈ H01 (U ), where ut is defined by the equation

(−∆ + 1)ut = et χ0 − χt (Equality in H −1 (U )),

(11a)

satisfies the following conditions, ut ≥ 0,

(11b)

hut , 1 − χt i = 0

(11c)

for every t ∈ [0, T ]. To motivate this definition, we prove (formally) that if {Bt }0≤t ≤T is a classical solution of (P), then {χt }0≤t ≤T is a weak solution of (P). The first step is to extend the definition of the pressure to all of U: p¯ t (x) = pt (x) for x ∈ Bt , and = 0 otherwise. Since pt ∈ C 1 (Bt ) (because Γt is C 2 ) and pt = 0 on Γt a standard result about Sobolev spaces tells us that p¯ t ∈ H01 (U ), Brézis [2, Proposition IX.18]. For any test function φ ∈ C0∞ (Rn ), integration by parts gives

Z h−∆p¯ t , φi = h¯pt , −∆φi = − pt ∆φ dx Bt Z Z Z ∂ pt = ∇ pt · ∇φ dx = φ dst − φ ∆pt dx, ∂n Bt Γt Bt where dst is the surface measure on Γt . It now follows from (1) and (2) that p¯ t satisfies the equation (−∆ + 1)¯pt = χt − Γ˙ t dst ,

(12)

which holds in the sense of distributions. To proceed further we need to express the normal velocity in terms of the characteristic functions χt . Recall that there exists a C 2 -function f = f (x, t ) such that Bt = {x ∈ U : f (x, t ) > 0}, so if H = H (s) denotes the Heaviside function, then

χt (x) = H (f (x, t )).

(13)

N.Chr. Overgaard / Nonlinear Analysis 70 (2009) 3658–3664

3661

If we take gradients on both sides of this identity, we get −n dst = ∇χt = H 0 (f )∇x f = δ0 (f )∇x f , where n is the outward unit normal of Bt and δ0 (f ) the pullback by f of the Dirac measure δ0 on R (well-defined because zero is a regular value for f ). Taking norms on both sides of this identity gives: dst = δ0 (f )|∇x f |.

(14)

The time derivative of (13), combined with (6) and (14), yields d dt

χt = δ 0 ( f )

∂f ∂ f /∂ t = |∇x f |δ0 (f ) = Γ˙ t dst . ∂t |∇x f |

(15)

Inserting this identity into (12) gives the pivotal equation of this section:

(−∆ + 1)¯pt = χt −

d dt

χt .

(16)

Now, notice that es−t is an integrating factor for the right-hand side of (16). If we introduce a new dependent variable Rt ut = ut (x) defined by the weighted Baiocchi transform ut = 0 et −s p¯ s ds (also used in [8], as it turns out), we find:

(−∆ + 1)ut =

Z

t

et −s (−∆ + 1)¯ps ds

0

Z =

t

et − s



χs −

0

d ds

χs



t

Z ds =

d ds

0


−et −s χs ds = et χ0 − χt ,

hence {χt }0≤t ≤T satisfies (11a). Since p¯ t ≥ 0 for all t, (11b) holds too, Finally, since the sets Bt are monotonically increasing, it is easy to see that (11c) is also satisfied. Thus {χt }0≤t ≤T is a weak solution of (P), as required. 5. Reformulation as a variational inequality and existence One cannot simply solve for ut in (11a) and set χt = χ{ut >0} , to obtain a weak solution to (P), because the unknown χt appears on the right-hand side of the equation. However, by essentially throwing away the χt -term, we can show that ut must be the solution of a variational inequality in the Hilbert space H01 (U ). Existence and uniqueness for this problem is well-known. In fact, it can be based almost entirely on the projection theorem. More precisely, suppose {χt }0≤t ≤T is a weak solution of (P). Fix t ∈ [0, T ], set ρt = et χ0 − 1, and rewrite (11a) in the form

(−∆ + 1)ut = ρt + (1 − χt ).

(17)

If we drop 1 − χt ≥ 0 and use (11b) and (11c), we find that ut is a solution of the following linear complementary problem:

(−∆ + 1)ut ≥ ρt , ut ≥ 0, h(−∆ + 1)ut − ρt , ut i = 0.

(18a) (18b) (18c)

Let K = {v ∈ H01 (U ) : v ≥ 0 a.e. on U }, then K is non-empty, closed and convex. It is well-known that ut solves (18a)–(18c) if and only if it solves the variational inequality and (ut , v − ut )H 1 ≥ hρt , v − ut i

ut ∈ K

0

where (u, v)H 1 = 0

R U

for all v ∈ K ,

(∇ u · ∇v + uv) dx is the scalar product on

H01

(19)

(U ) and h·, ·i the duality bracket for

H01

(U ) and H

−1

(U ).

This reformulation of the original problem is the key to the existence and uniqueness for solutions to (P): Proposition 4. For each t ≥ 0 the variational inequality (19) has a unique solution ut ∈ H01 (U ). The mapping [0, ∞) 3 t 7→ ut ∈ H01 (U ) is continuous and satisfies the bound kut kH 1 ≤ (et − 1)kχ0 kL2 . In particular u0 = 0. 0

Proof. We briefly recall the argument: By Riesz’ lemma there exists a unique vector ρˆ t ∈ H01 (U ) such that (ρˆ t , v)H 1 = hρt , vi for all v ∈ H01 (U ), (i.e. (−∆ + 1)ρˆ t = ρt in H −1 (U ).) Problem (19) thus becomes ut ∈ K

and (ρˆ t − ut , v − ut )H 1 ≤ 0 0

0

for all v ∈ K ,

or, equivalently, to find ut ∈ K such that kρˆ t − ut kH 1 = infv∈K kρˆ t − vkH 1 . The unique solution of this problem is ut = PK ρˆ t , 0 0 where PK denotes the projection onto K (see, e.g., Brézis [2, Théorème V.2]). Since PK does not increase the norm, we get continuity in time: For any s, t ∈ [0, T ],

kut − us kH 1 ≤ kρˆ t − ρˆ s kH 1 = kρt − ρs kH −1 ≤ kχ0 kL2 |et − es |. 0

0

The bound on kut kH 1 is obtained by setting s = 0 and using u0 = 0. The latter follows from (19) because ρ0 ≤ 0. 0



3662

N.Chr. Overgaard / Nonlinear Analysis 70 (2009) 3658–3664

It is a standard result that if ut solves (19) then ut ∈ W 2,p (U ) for all 1 ≤ p < ∞. Hence ut ∈ C 1,α (U ), 0 < α < 1, by Sobolev’s inequality, see [9] or the simplified proof in Gustafsson [7]. From now on, a map {ut }0≤t ≤T is called a variational solution of (P) with initial data χ0 , if ut solves (19) with ρt = et χ0 − 1, and has compact support in U, for all t ∈ [0, T ]. To analyze this concept we need the following comparison result, see [6, Lemma 4.1]. Proposition 5. Let ut be the solution of the variational inequality (19). If v ∈ K satisfies (−∆ + 1)v ≥ ρt , then ut (x) ≤ v(x) a.e. in U. Let u denote the solution of (19) when ρt is replaced by ρ ∈ H −1 (U ). If ρ ≥ ρt , as distributions, then the proposition implies that u ≥ ut a.e. in U. This principle is used in the proof of the following result, which shows that variational solutions of (P) always exist. (Here U is a ball.) Corollary 6. (a) If {ut }0≤t <∞ is the map in Proposition 4, then ut ≥ us for 0 ≤ s ≤ t < ∞. (b) Assume supp χ0 is contained in the ball B0 = {x : |x| < R0 }. Let {BRt }0≤t ≤T be the radially symmetric classical solution of (P) with BR0 = B0 , where T > 0 is such that the closure of BRt is contained in U, for all t ∈ [0, T ]. Then supp ut ⊂ BRt for all t ∈ [0, T ]. In particular, {ut }0≤t ≤T is a variational solution of (P). Proof. (a) follows from the observation ρt ≥ ρs for 0 ≤ s ≤ t < ∞. To prove (b), set ρtR = et χBR − 1 and let uRt denote the t

corresponding solution of (19). By assumption, ρt ≤ ρtR hence ut ≤ uRt for all t ∈ [0, T ]. The assertion now follows because supp uRt = BRt , by Theorem 2 and the construction of weak solutions from classical ones in Section 4.  6. From variational to weak solutions It is now shown that every variational solution of (P) is also a weak solution of (P), that is, the two notions of solvability are equivalent. Theorem 7. Suppose {ut }≤t ≤T is a variational solution of (P) with initial data χ0 . For 0 < t ≤ T , define χt = χ{x:ut (x)>0} , and set χt = χ0 when t = 0. Then {χt }0≤t ≤T is a weak solution of (P). Proof. Since each ut is continuous (regularity) with compact support (by definition), Bt = {ut > 0} is open and bounded, so the map [0, T ] 3 t 7→ χt ∈ L2 (U ) is well-defined. It remains to check the conditions (11a)–(11c). First, observe that since Bt is open and ut ∈ H 2 (U ), then (18a) and (18c) imply

(−∆ + 1)ut =

 ρt 0

a.e. on Bt a.e. on U \ Bt = {ut = 0}.

Therefore, (−∆ + 1)ut = χt ρt = χt (et χ0 − 1) = et χ0 − χt on U because χt ≥ χ0 , by part (a) of Corollary 6). Thus (11a) holds. By the above computation, (−∆ + 1)ut − ρt = 1 − χt , so (18c) implies (11c). Finally, (18b) is identical to (11b), so the proof is complete.  Since the choice the encapsulating set U, used in the construction of variational solutions, is inessential, Theorem 7 implies the following existence result: Corollary 8. To each bounded open set B0 in Rn there exists a unique mapping {χt }0≤t <∞ such that: (a) each χt is the characteristic function of a bounded open set and χ0 = χB0 ; (b) t 7→ χt is increasing; (c) For any T > 0, the restricted map {χt }0≤t ≤T is a weak solution of (P) with initial data χB0 ; (d) If B0 ⊂ {x : |x − a| < R0 }, for some a ∈ Rn and R0 > 0, then χt ≤ χ{x:|x−a| 0. In particular, Bt will not undergo changes in topology. Such a result is true for the Hele-Shaw flow with injection through a star-shaped core, as proved by Di Benedetto and Friedman [1]. However the proof given there cannot be immediately modified to our problem, which remains open. Acknowledgements The author wishes to thank Björn Gustafsson at KTH, Stockholm, and Professor Anders Melin at Lund University for helpful comments and discussions. Appendix. Proof of Proposition 1 To analyze Gν (z ) = (iz /2)−ν Jν (iz ), r ∈ R, defined in (8), recall Lommel’s integral representation of the Bessel function Jν (z ), cf. Watson [12, p. 38]:

N.Chr. Overgaard / Nonlinear Analysis 70 (2009) 3658–3664

3663

Fig. 1. A solution of (P) obtained by solving the linear complementary problem (18a)–(18c) using the PSOR algorithm. The initial biofilm colony B0 is the star-shaped region shown in black, and the surrounding closed curves correspond to the biofilm interfaces Γt at time t = k∆t, k = 1, . . . , 12, and time step ∆t = 0.25.

( 12 z )ν

Jν (z ) =

π

Z

Γ (ν + 21 )Γ ( 21 )

cos(z cos θ ) sin2ν θ dθ .

0

If we make the change of variables t = − cos θ , employ symmetry, and finally set z = ir, then we find: Gν (r ) =

Z

1

Γ (ν + )Γ ( ) 1 2

1 2

1

1

(1 − t 2 )ν− 2 ert dt .

−1

This integral is convergent if Re ν > L -function f (t ) = (1 − 1

ν− 1 t 2 )+ 2 .

1 , 2

and is real-valued if ν ∈ R. We see that Gν is essentially the Laplace transform of the

Therefore Proposition 1 is a consequence of the result:

Theorem 9. Suppose f ∈ L1 (R) is non-negative almost everywhere with compact support (in distributional sense), and let F (r ) be defined by the integral F (r ) =

Z

∞

f (s)ers ds,

(r ∈ R).

(A.1)

−∞

If f is not identically zero in L1 (R), then the function F 0 (r )/F (r ) is strictly increasing on R. (In particular, F (r ) is strictly logarithmically convex.) Moreover, if [a, b] is the smallest interval which contains supp f , then lim F 0 (r )/F (r ) = b

r →∞

and

lim F 0 (r )/F (r ) = a.

(A.2)

r →−∞

rt Proof. R ∞Define w(r ; t ) = f (t )e /F (r ). For each r ∈ R, w(r ; ·) is a weight function because w(r ; t ) R≥∞0, supp w(r ; ·) ⊂ [a, b], and −∞ w(r ; t ) dt = 1. Differentiation under the integral sign in (A.1) gives F 0 (r )/F (r ) = −∞ t w(r ; t ) dt, so clearly a < F 0 (r )/F (r ) < b. Differentiating once more yields

d F 0 (r ) dr F (r )

Z

∞

t w(r ; t ) dt − 2

=

Z

−∞

∞

t w(r ; t ) dt

2

> 0,

−∞

by the Schwarz inequality. Equality in this estimate would imply that t = 1 on supp w(r ; t ), which is impossible. Thus F 0 (r )/F (r ) is strictly increasing on R, and it only remains R ∞ to determine its limits at infinity. In order to compute limr →∞ F 0 (r )/F (r ), observe that the definition of b implies b−/2 f (s) ds > 0, for any  > 0, hence

Z

R b−

b−

w(r ; t ) dt = −∞

R−∞ ∞ −∞

er (b−) −∞ f (t )dt −∞ f (t )dt R∞ ≤ ≤ e−r /2 R ∞ . rs r ( b −/ 2 ) f (s)e ds e f (s)ds f (s)ds b−/2 b−/2 f (t )ert dt

R b−

R b−

This shows that the mass of w(r ; ·) on the interval (−∞, b − ] tends to zero rapidly as r → ∞. For any  > 0, and r sufficiently large, it now follows that

3664

N.Chr. Overgaard / Nonlinear Analysis 70 (2009) 3658–3664

0
F 0 (r ) F (r )

Z

∞

(b − t )w(r ; t ) dt =

=

b−

Z

Z

b−

w(r ; t ) dt +  −∞

∞

(b − t )w(r ; t ) dt

+ b−

−∞

−∞

≤ (b − a)

Z

Z

∞

w(r ; t ) dt ≤ , b−

proving the first limit in (A.2). The second limit can be verified similarly.



G0ν (r )/Gν (r )

Since = Jν+1 (ir )/iJν (ir ), the limits in Proposition 1 follow from the asymptotic formulas for the Bessel functions [12, Chapter VII]. However, this does not yield the strict monotonicity of G0ν (r )/Gν (r ). References [1] E.D. Benedetto, A. Friedman, The ill-posed Hele-Shaw model and the Stefan problem for supercooled water, Trans. Amer. Math. Soc. 282 (1984) 183–204. [2] H. Brézis, Analyse fonctionelle-Theorie et applications, Dunod, Paris, 1999. [3] J. Crank, Free and Moving Boundary Problems, Oxford University Press, 1984. [4] J. Dockery, I. Klapper, Finger formation in biofilm layers, SIAM J. Appl. Math. 62 (3) (2001) 853–869. [5] C.M. Elliott, V. Janovský, A variational inequality approach to Hele-Shaw flow with a moving boundary, Proc. Roy. Soc. Edinburgh 88A (1981) 93–107. [6] B. Gustafsson, Application of variational inequalities to a moving boundary problem for Hele Shaw flows, SIAM J. Math. Anal. 16 (1985) 279–300. [7] B. Gustafsson, A simple proof of the regularity theorem for the variational inequality of the obstacle problem, Nonlinear Anal. TMA 10 (1986) 1487–1490. [8] B. Gustafsson, Direct and inverse balayage—Some new developments in classical potential theory, Nonlinear Anal. TMA 30 (1997) 2557–2565. [9] D. Kinderlehrer, G. Stampacchia, An Introduction to Variational Inequalities and their Applications, in: Classics in Applied Mathematics, SIAM, 2000. [10] J.C. Kissel, P.L. McCarty, R.L. Street, Numerical simulation of mixed-culture biofilm, J. Environmental Eng. 110 (1984) 393–411. [11] J.R. Ockendon, S.D. Howison, A.A. Lacey, Mushy regions in negative squeeze films, Q. J. Mech. Appl. Math. 56 (2003) 361–379. [12] G.N. Watson, A Treatise on the Theory of Bessel Functions, second ed., Cambridge University Press, 1996.