Skip to content
BY 4.0 license Open Access Published by De Gruyter May 27, 2020

Weak and stationary solutions to a Cahn–Hilliard–Brinkman model with singular potentials and source terms

Matthias Ebenbeck and Kei Fong Lam


We study a phase field model proposed recently in the context of tumour growth. The model couples a Cahn–Hilliard–Brinkman (CHB) system with an elliptic reaction-diffusion equation for a nutrient. The fluid velocity, governed by the Brinkman law, is not solenoidal, as its divergence is a function of the nutrient and the phase field variable, i.e., solution-dependent, and frictionless boundary conditions are prescribed for the velocity to avoid imposing unrealistic constraints on the divergence relation. In this paper we give a first result on the existence of weak and stationary solutions to the CHB model for tumour growth with singular potentials, specifically the double obstacle potential and the logarithmic potential, which ensures that the phase field variable stays in the physically relevant interval. New difficulties arise from the interplay between the singular potentials and the solution-dependent source terms, but can be overcome with several key estimates for the approximations of the singular potentials, which maybe of independent interest. As a consequence, included in our analysis is an existence result for a Darcy variant, and our work serves to generalise recent results on weak and stationary solutions to the Cahn–Hilliard inpainting model with singular potentials.

MSC 2010: 35K35; 35D30; 35J61; 35Q92; 92C50; 76D07

1 Introduction

Phase field models [17, 45] have recently emerged as a new mathematical tool for tumour growth, which offer new advantages over classical models [11, 24] based on a free boundary description, such as the ability to capture metastasis and other morphological instabilities like fingering in a natural way. In conjunction with clinical data, statistical methodologies and image analysis, they have begun to impact the development of personalised cancer treatments [5, 6, 37, 39, 41, 42].

To support this continuing effort, in this paper we expand the recent analysis performed in [19, 20] for a class of phase field tumour models based on the Cahn–Hilliard–Brinkman (CHB) system. For a bounded domain Ω ⊂ ℝd for d ∈ {2, 3} with boundary Σ : = Ω and outer unit normal n, and an arbitrary but fixed terminal time T, we consider for ΩT := Ω × (0, T) the following system of equations

div(v)=Γv(φ,σ):=bv(φ)σ+fv(φ) in ΩT, (1.1a)
div(T(v,p))+νv=(μ+χσ)φ in ΩT, (1.1b)
tφ+div(φv)Δμ=Γφ(φ,σ):=bφ(φ)σ+fφ(φ) in ΩT, (1.1c)
μ=ψ(φ)Δφχσ in ΩT, (1.1d)
0=Δσh(φ)σ in ΩT, (1.1e)

where the viscous stress tensor T and the symmetrised velocity gradient Dv are given by


The model (1.1) is a description of the evolution of a two-phase cell mixture, containing tumour cells and healthy host cells, surrounded by a chemical species acting as nutrients only for the tumour cells, and is transported by a fluid velocity field. The variable σ denotes the concentration of the nutrient and is assumed to evolve quasi-statically, while φ denotes the difference in the volume fractions of the cells, with the region {φ = 1} representing the tumour cells and {φ = −1} representing the host cells. The fluid velocity v is taken as the volume averaged velocity, with pressure p, while μ denotes the chemical potential associated to φ.

Equations (1.1c)-(1.1d) comprise a convective Cahn–Hilliard system for (φ, μ) with source term Γφ(φ, σ), that couples the nutrient equation (1.1e) also through the consumption term h(φ) σ. The constant parameter χ appearing in (1.1d) captures the chemotaxis effect, see [32]. On the other hand, equations (1.1a)-(1.1b) form a Brinkman system for (v, p) with the term (μ + χ σ) ∇ φ modelling capillary forces, and in the definition of the stress tensor, the functions η and λ represent the shear and bulk viscosities, respectively. An interesting contrast with previous phase field models in two-phase flows that employ a volume averaged velocity, such as [4, 10, 18], is that the velocity in (1.1) is not solenoidal. This can be attributed to the fact that in the case of unmatched densities, the gain σ bφ(φ) and loss fφ(φ) of cellular volume leads to sources σ bv(φ) and sinks fv(φ) in the mass balance, see for example [32, 40] for more details in the model derivation.

As a consequence of (1.1a), we find the relation


and the typical no-penetration boundary condition vn = 0 on Σ × (0, T) will lead to a mean-zero compatibility condition for Γv(φ, σ), which can not be fulfilled in general since Γv depends on φ and σ. For models with Darcy’s law instead of (1.1b), there are some works in the literature [25, 30] prescribing alternative boundary conditions, such as zero transport flux ∇ μn = φ vn and zero pressure p = 0 on ΣT := Σ × (0, T) to circumvent the compatibility condition on Γv. For the Brinkman model, the compatibility issue can be avoided by prescribing the frictionless boundary condition T n = 0 on ΣT, as considered in [19, 20, 21, 22].

Hence, we furnish (1.1) with the following initial and boundary conditions

nμ=nφ=0 on ΣT, (1.2a)
nσ=K(1σ) on ΣT, (1.2b)
T(v,p)n=0 on ΣT, (1.2c)
φ(0)=φ0 in Ω, (1.2d)

where n f = ∇ fn is the normal derivative on Σ, φ0 is a prescribed initial datum and K is a positive boundary permeability constant.

The resulting system (1.1)-(1.2) has been studied previously by the first author in a series of works [19, 20, 21, 22] concerning well-posedness, asymptotic limits and optimal control, in the setting where the double well potential ψ, whose first derivative appears in (1.1d), is continuously differentiable. The common example is the quartic potential ψ(s) = (s2 − 1)2. The purpose of the present paper is to provide a first result for the weak existence to the CHB model when ψ is a singular potential, where either the classical derivative of ψ does not exist, or the derivative ψ′ blows up outside a certain interval. For the former, the prime example is the double obstacle potential [8, 9]

ψdo(r)=12(1r2)+I[1,1](r)=12(1r2) for r[1,1],+ otherwise,

and for the latter the logarithmic potential

ψlog(r)=θ2((1+r)log(1+r)+(1r)log(1r))+θc2(1r2) for r(1,1)

with positive constants 0 < θ < θc is the standard example.

Let us briefly motivate the need to consider singular potentials such as ψdo and ψlog. The variable φ has a physical interpretation as the difference between the volume fractions of the tumour cells and the host cells. As volume fractions are non-negative and bounded above by 1, φ should lie in the physical range [−1, 1]. This natural boundedness of φ becomes particularly important for physical models in fluid flow and biological sciences, where the mass density of the mixture expressed as an affine linear function of φ may become negative if φ strays out of [−1, 1]. A counterexample in [14] shows that for polynomial ψ, the phase field variable φ can take values outside [−1, 1], and hence a remedy to enforce the natural bounds for φ is to introduce non-smooth potentials into the model.

For the original Cahn–Hilliard equation and its variants in solenoidal two-phase flow with the boundary condition (1.2a), there are numerous contributions in the literature involving singular potentials, of which we cite [1, 3, 8, 16, 27, 35, 36] and refer to the references cited therein. The most challenging aspect with singular potentials in the global existence of weak solutions is to show that the chemical potential μ is controlled in the Bochner space L2(0, T; H1(Ω)). This is equivalent to controlling the mean value μΩ := 1|Ω| Ω μ dx in L2(0, T), since the gradient ∇ μ is controlled from basic energy identities. For polynomial ψ this can be done via the relation Ω μ dx = Ω ψ′(u) dx with suitable growth conditions on ψ′, but this argument fails for the singular case. The techniques used in the aforementioned references depend on first establishing the assertion φΩ(t) ∈ (−1, 1) for all t > 0, which holds automatically for suitably chosen initial data thanks to the property of mass conservation φΩ(t) = φΩ(0) for all t > 0.

When the Cahn–Hilliard component is affixed with a source term, like (1.1c), then in general the mass conservation property is lost. Therefore, simply using previous techniques would only give a local-in-time weak existence result holding as long as φΩ(t) ∈ (−1, 1), see for instance the work of [13] on the so-called Cahn–Hilliard inpainting model [7] with logarithmic potential obtained from (1.1c)-(1.1d) by setting ψ = ψlog, v = 0, σ = 0 and fφ(φ) = λ (Iφ) for given non-negative λL(Ω) and ∣I∣ ≤ 1 a.e. in Ω. Recently, the second author has established a global existence result to the inpainting model with the double obstacle potential ψdo in [33]. The key ingredients involved are a careful analysis of the interplay between the source terms and approximations of the singular part 𝕀[−1,1] of ψdo, as well as showing that the source terms lead to the conclusion that φΩ(t) ∈ (−1, 1) for all t > 0. Then, previous methods for controlling μΩ can be applied to achieve global weak existence.

To generalise the methodology of [33] for the CHB model (1.1) with ψdo, we found that it suffices to prescribe suitable conditions on bv, bφ, fv and fφ at φ = ± 1. Namely, we ask that bv and bφ are non-negative and vanish at ± 1, while fφfv is negative (resp. positive) at 1 (resp. −1). In Remark 2.1 below we give a biologically relevant example of bv, bφ, fv and fφ satisfying the above conditions. With the help of a modified version of a key estimate (see Proposition 3.3), the same methodology can be used to show global weak existence for the CHB model (1.1) with ψlog, thereby improving also the results of [13, 43] concerning the inpainting model with the logarithmic potential ψlog. Uniqueness for the CHB tumour model remains an open problem, as the frictionless boundary condition (1.2c) and the presence of source terms Γv and Γφ introduce new difficulties that do not permit us to mimic the arguments developed in [16, 34, 35].

By formally sending the bulk and shear viscosities in the stress tensor T to zero, the CHB model (1.1) reduces to a Cahn–Hilliard–Darcy (CHD) tumour model with


replacing (1.1b). Without source terms and the nutrient, the CHD model is also called the Hele-Shaw–Cahn–Hilliard system, where recent progress on strong well-posedness with the logarithmic potential can be found in [34, 35]. Meanwhile, for the CHD tumour model with polynomial ψ, it appears that the regularity of the weak solutions in [29, 30] is insufficient to replicate the usual procedure of approximating the singular potential with a sequence of polynomial potentials, deriving uniform estimates and passing to the limit. To the authors’ best knowledge the global weak existence to the CHD tumour model with singular potentials remains an open problem, and inspired by the relationship between the Brinkman and Darcy laws, we employ an idea of [20] to deduce the global weak existence to the CHD tumour model with singular potentials by scaling the shear and bulk viscosities in an appropriate way.

The second contribution of the present paper is to provide a first result for stationary solutions to the CHB model with singular potentials, where the time derivative in (1.1c) vanishes but the convection term remains, i.e., we allow for the possibility of a non-zero velocity field. This is interesting as the fluid velocity can be used to inhibit the growth of the tumour cells, and in an optimal control framework, these stationary solutions are ideal candidates in a tracking-type objective functional, e.g. see [12]. Unfortunately, the particular form of the source terms Γv(φ, σ) and Γφ(φ, σ) indicates that the CHB model does not admit an obvious Lyapunov structure, and so it is unlikely that the a priori estimates used to prove weak existence are uniform in time. This is in contrast to the phase field tumour model of [38] studied in [12, 15, 26, 48, 49] where tools such as global attractors and the Łojasiewicz–Simon inequality can be used to quantify the long-time behaviour of global solutions. We also mention the recent work [44] establishing a global attractor for a simplified model similar to the subsystem (1.1c)-(1.1e) with v = 0 and χ = 0, where despite the lack of an obvious Lyapunov structure, the authors can derive dissipative estimates under some constraints on the model parameters.

However, we opt to prove directly the existence of stationary solutions using the methodology of [33], as opposed to investigate the long-time behaviour of time-dependent solutions. Thanks to the well-posedness of the Brinkman subsystem (1.1a), (1.1b), (1.2c) and the nutrient subsystem (1.1e), (1.2b), we can express σ = σφ, v = vφ and p = pφ as functions of φ, since μ can also be interpreted as a function of φ by (1.1d). Then, the stationary CHB system is formally equivalent to the fourth order elliptic problem

div(φvφ)+Δ2φΔψ(φ)+χΔσφ=Γφ(φ,σφ) in Ω,nφ=n(Δφ+χσφ)=0 on Σ.

The novelty of our work lies in controlling the convection term div(φ vφ) with the help of the unique solvability of the Brinkman system, and the existence result (Theorem 2) is a non-trivial extension of [33]. Furthermore, this also shows the existence of stationary solutions to the inpainting model with logarithmic potential, thereby completing the analysis of weak and stationary solutions to the Cahn–Hilliard inpainting model with both singular potentials. We mention that a similar strategy does not work for the stationary Cahn–Hilliard–Darcy system with non-zero velocity field, as regularity of the corresponding velocity field is simply insufficient.

The remainder of this paper is organised as follows. In section 2 we cover several useful preliminary results and state our main results on the existence of global weak solutions and stationary solutions to the CHB model (1.1)-(1.2) with singular potentials, as well as the global weak existence to the CHD tumour model. The proofs of these results are then contained in sections 3, 4 and 5, respectively. Proofs of several useful properties of approximations to both singular potentials are contained in appendix A, while in appendix B we give a proof of a well-posedness result for the Brinkman system that plays a significant role in our work.

2 Main results

2.1 Notation and preliminaries

For a real Banach space X we denote by ∥⋅∥X its norm, by X* its dual space and by 〈⋅, ⋅〉X the duality pairing between X* and X. If X is an inner product space the associated inner product is denoted by (⋅, ⋅)X. For two matrices A = (ajk)1≤j,kd, B = (bjk)1≤j,kd ∈ ℝd×d, the scalar product A : B is defined as the sum j,k=1d ajkbjk.

For 1 ≤ p ≤ ∞, r ∈ (1, ∞), β ∈ (0, 1) and k ∈ ℕ, the standard Lebesgue and Sobolev spaces on Ω and on Σ are denoted by Lp := Lp(Ω), Lp(Σ), Wk,p := Wk,p(Ω), and Wβ,r(Σ) with norms ∥⋅∥Lp, ∥⋅∥Lp(Σ), ∥⋅∥Wk,p and ∥⋅∥Wβ,r(Σ), respectively. In the case p = 2 we use the notation Hk := Wk,2 with the norm ∥⋅∥Hk. The space H01 is the completion of C0 (Ω) with respect to the H1-norm, while the space Hn2 is defined as {wH2 : n w = 0 on Σ}. For the Bochner spaces, we use the notation Lp(X) := Lp(0, T; X) for a Banach space X with p ∈ [1, ∞], and


for two or more Bochner spaces A and B. For the dual space X* of a Banach space X, we introduce the (generalised) mean value by

vΩ:=1|Ω|Ωvdxfor vL1,vΩ:=1|Ω|v,1Xfor vX,

and also the subspace of L2-functions with zero mean value:


The corresponding function spaces for vector-valued or tensor-valued functions are denoted in boldface, i.e., Lp, Wk,p, Hk, Lp(Σ) and Wβ,r(Σ). Furthermore, we introduce the space


equipped with the norm


where div is the weak divergence. We now state three auxiliary lemmas that are used for the study of our model. The first lemma concerns the solvability of the divergence problem, which will be useful for the mathematical treatment of equation (1.1a).

Lemma 2.1

([28, Sec. III.3]). Let Ω ⊂ ℝd, d ≥ 2, be a bounded domain with Lipschitz-boundary and let 1 < q < ∞. Then, for every fLq and aW1−1/q,q(Σ) satisfying

Ωfdx=ΣandH, (2.1)

there exists at least one solution uW1,q to the problem

div(u)=finΩ,u=aonΣ, (2.2)

satisfying for some positive constant C = C(Ω, q) the estimate

uW1,qCfLq+aW11/q,q(Σ). (2.3)

Since we invoke Lemma 2.1 frequently with a=1|Σ|(Ωfdx)n, we introduce

u=D(f)div(u)=f in Ω,u=1|Σ|(Ωfdx)n on Σ,

where by the smoothness of the normal n, see (A1) below, (2.3) yields


The second lemma concerns the existence of weak solutions to the model (1.1) which forms the basis of our approximation procedure.

Lemma 2.2

([20, Thm. 2.5]). Suppose the following assumptions are satisfied:

  1. For d ∈ {2, 3}, Ω ⊂ ℝd is a bounded domain with C3-boundary.

  2. The positive constants ν, K and the non-negative constant χ are fixed.

  3. The viscosities η and λ belong to C2(ℝ) ∩ W1,∞(ℝ) and satisfy

    η0η(t)η1,0λ(t)λ0tR, (2.4)

    for positive constants η0, η1 and a non-negative constant λ0. The function hC0(ℝ) ∩ L(ℝ) is non-negative

  4. The source terms Γv and Γφ are of the form

    Γv(φ,σ)=bv(φ)σ+fv(φ),Γφ(φ,σ)=bφ(φ)σ+fφ(φ), (2.5)

    where bv, fvC0,1(ℝ) ∩ W1,∞(ℝ) and bφ, fφC0(ℝ) ∩ L(ℝ).

  5. The initial condition φ0 belongs to H1.

  6. The function ψC2(ℝ) is non-negative and satisfies

    ψ(s)R0|s|2R1,|ψ(s)|R21+|s|,|ψ(s)|R3sR, (2.6)

    where R0, R1, R2, R3 are positive constants.

Then, there exists a quintuple (φ, μ, σ, v, p) satisfying


and is a weak solution to (1.1)-(1.2) in the sense that (1.1a) holds a.e. in ΩT and

0=Ω2η(φ)Dv:DΦ+(λ(φ)div(v)p)div(Φ)+νvΦ(μ+χσ)φΦdx, (2.7a)
0=tφ,ζH1+Ωμζ+(φv+φΓv(φ,σ)Γφ(φ,σ))ζdx, (2.7b)
0=Ω(μ+χσ)ζψ(φ)ζφζdx, (2.7c)
0=Ωσζ+h(φ)σζdx+ΣK(σ1)ζdH, (2.7d)

for a.e. t ∈ (0, T) and for all ΦH1 and ζH1. Furthermore, there exists a positive constant C not depending on (φ, μ, σ, v, p) such that

φH1((H1))L(H1)L4(H2)L2(H3)+σL(H2)+μL4(L2)L2(H1)+div(φv)L2(L2)+vL8/3(H1)+pL2(L2)C. (2.8)

Remark 2.1

We point out that the assumptions on the potential in [20] are more restrictive in the case of potentials with quadratic growth. However, it can be checked easily that the proof of [20, Thm. 2.5] follows through when using potentials that satisfy (2.6). Furthermore, in checking the proof of [20, Thm. 2.5], it turns out that Lipschitz continuity of bv and fv are sufficient for the assertions of Lemma 2.2, as opposed to the C1-regularity stated in [20, (A4)].

The third lemma concerns the solvability of the Brinkman system, which we will use to study stationary solutions. A proof is contained in appendix, and we mention that similar results can be found in [2] in the context of resolvent operators and H-calculus.

Lemma 2.3

Let Ω ⊂ ℝd, d = 2, 3, be a bounded domain with C3-boundary Σ and outer unit normal n. Let cW1,r with r > d be given and fix exponent qr such that q > 1 for d = 2 and q 65 for d = 3, and suppose the functions η(⋅) and λ(⋅) satisfy (A3). Then, for any fLq, gW1,q, hW1−1/q,q(Σ), there exists a unique solution (v, p) ∈ W2,q × W1,q to

div(2η(c)Dv+λ(c)div(v)I)+νv+p=fa.e.inΩ, (2.9a)
div(v)=ga.e.inΩ, (2.9b)
(2η(c)Dv+λ(c)div(v)IpI)n=ha.e.onΣ, (2.9c)

satisfying the following estimate

vW2,q+pW1,qC(fLq+gW1,q+hW11q,q(Σ)), (2.10)

with a constant C depending only on η0, η1, λ0, ν, q, ∥cW1,r and Ω.

2.2 The time-dependent problem

We begin with a suitable notion of weak solutions for the model with the double obstacle potential ψdo and the logarithmic potential ψlog.

Definition 2.1

A quintuple (φ, μ, σ, v, p) is a weak solution to the CHB system (1.1)-(1.2) with the double obstacle potential ψdo if the following properties hold:

  1. the functions satisfy


    with φ(0) = φ0 a.e. in Ω.

  2. equation (1.1a) holds a.e. in ΩT, while (2.7a), (2.7b) and (2.7d) are satisfied for a.e. t ∈ (0, T) and for all ΦH1 and ζH1.

  3. for a.e. t ∈ (0, T), φ(t) ∈ 𝓚 := { fH1 : ∣f∣ ≤ 1 a.e. in Ω} and

    Ω(μ+χσ+φ)(ζφ)φ(ζφ)dx0ζK. (2.11)

We say that (φ, μ, σ, v, p) is a weak solution to (1.1)-(1.2) with the logarithmic potential ψlog if properties (a) and (b) hold along with

  1. φ(x, t)∣ < 1 for a.e. (x, t) ∈ ΩT and for a.e. t ∈ (0, T),

    μ+χσψlog(φ)+Δφ=0a.e.inΩT. (2.12)

Our first result concerns the existence of weak solutions to the CHB system (1.1)-(1.2) with singular potentials.

Theorem 1

Suppose (A1)-(A3) hold along with

  1. The source terms Γv and Γφ are of the form (2.5) with fvC0,1([−1, 1]), fφC0([−1, 1]), bvC0,1([−1, 1]), bφC0([−1, 1]) satisfying

    bv(±1)=bφ(±1)=0,fφ(1)fv(1)<0,fφ(1)+fv(1)>0. (2.13)
  2. The initial condition φ0 belongs to 𝓚 with ∣(φ0)Ω∣ < 1.

Then, there exists a weak solution (φ, μ, σ, v, p) to (1.1)-(1.2) with the double obstacle potential ψdo in the sense of Definition 2.1, and additionally satisfies

0σ1a.e.inΩT,φL2(0,T;W2,p)L4(0,T;Hn2) (2.14)

for any p ∈ [2, ∞) if d = 2 and p = 6 if d = 3.

In addition, if the following assumptions are satisfied:

  1. There exists a constant c such that for any 0 < δ ≪ 1,

  2. bφ(s) log(1+s1s) and bv(s) log(1+s1s) C0([−1, 1]).

  3. The initial condition φ0H1(Ω) satisfiesφ0(x)∣ < 1 for a.e. xΩ and β̂log(φ0) ∈ L1(Ω).

Then, there exists a weak solution (φ, μ, σ, v, p) to (1.1)-(1.2) with the logarithmic potential ψlog in the sense of Definition 2.1 and satisfies (2.14). Furthermore, for a.e. t ∈ (0, T), it holds that

Ωψlog(φ(t))+12|φ(t)|2dx+0tΩ|μ|2+2η(φ)|Dv|2+ν|v|2dxdt0tΩχμσ+(ΓφφΓv)(ψlog(φ)Δφ)dxdt+0tΩ2η(φ)Dv:Du+νvudx(ψlog(φ)Δφ)φudxdt+Ωψlog(φ0)+12|φ0|2dx, (2.15)

where u = 𝓓(Γv(φ, σ)).

Example 2.1

We now give a biologically relevant example of the source terms Γv and Γφ that satisfy (2.13). Following [32, Sec. 2.4.1], in a domain Ω occupied by both tumour cells and healthy cells, we denote by ρ1 the actual mass density of the healthy cells per unit volume in Ω and by ρ̄1 the (constant) mass density of the healthy cells occupying the whole of Ω. Then, it follows that ρ1 ∈ [0, ρ̄1] and the volume fraction of the healthy cells can be defined as the ratio u1=ρ1ρ¯1 ∈ [0, 1]. Let ρ2, ρ̄2 and u2 be the actual mass density of the tumour cells per unit volume in Ω, the (constant) mass density of the tumour cells occupying the whole of Ω, and the volume fraction of the tumour cells, respectively. Assuming there are no external volume compartment aside from the tumour and healthy cells, we have u1 + u2 = 1. Then, setting φ := u2u1, we consider the function


where P, A > 0 are the constant proliferation and apoptosis rates, so that proliferation occurs only at the interface region {−1 < φ < 1}. From the modelling, we have


and so we can set


where bv(± 1) = bφ(± 1) = 0 and


It is also clear that bφ(s) and bv(s) satisfy (C1) and (C2). Furthermore, we point out that it is possible to have α < 0 in the case ρ̄2 > ρ̄1.

2.3 The stationary problem

The stationary problem comprises of the equations (1.1a), (1.1b), (1.1d), (1.1e) posed in Ω and

div(φv)=Δμ+Γφ(φ,σ)in Ω, (2.16)

together with the boundary conditions (1.2a)-(1.2c) posed on Σ.

Definition 2.2

A quintuple (φ, μ, σ, v, p) is a weak solution to the stationary CHB system with the double obstacle potential ψdo if the following properties hold:

  1. the functions satisfy

  2. equation (1.1a) holds a.e. in Ω, while (2.7a), (2.7d) and

    0=Ωμζ+(φv+φΓv(φ,σ)Γφ(φ,σ))ζdx (2.17)

    are satisfied for all ΦH1 and ζH1.

  3. (2.11) holds along with φ ∈ 𝓚 = {fH1 : ∣f∣ ≤ 1 a.e. in Ω}.

We say that (φ, μ, σ, v, p) is a weak solution to the stationary CHB system with the logarithmic potential ψlog if properties (d) and (e) are satisfied and

  1. (2.12) holds along withφ(x)∣ < 1 for a.e. xΩ.

Theorem 2

Under (A1)-(A3) and (B1), there exists a weak solution (φ, μ, σ, v, p) to the stationary CHB model with double obstacle potential ψdo in the sense of Definition 2.2. If, in addition, (C1) and (C2) hold, then there exists a weak solution (φ, μ, σ, v, p) to the stationary CHB model with logarithmic potential ψlog in the sense of Definition 2.2. Furthermore, in both cases, 0 ≤ σ ≤ 1 a.e. in ΩT and the following regularity properties hold


for q ∈ [2, ∞) if d = 2 and q = 6 if d = 3. In particular, the quintuple (φ, μ, σ, v, p) is a strong stationary solution.

2.4 Darcy’s law

By setting the viscosities η(⋅) and λ(⋅) to zero, the CHB model (1.1)-(1.2) reduces to a Cahn–Hilliard–Darcy (CHD) model consisting of (1.1a), (1.1c)-(1.1e) and the Darcy law

v=1ν(p(μ+χσ)φ)in ΩT, (2.18)

furnished with the initial-boundary conditions (1.2a)-(1.2b), (1.2d) together with the Dirichlet boundary condition

p=0on ΣT. (2.19)

We begin with a notion of weak solutions for the CHD model with singular potentials.

Definition 2.3

A quintuple (φ, μ, σ, v, p) is a weak solution to the CHD system (1.1a), (1.1c)-(1.1e), (1.2a)-(1.2b), (1.2d), (2.18), (2.19) with the double obstacle potential ψdo if property (c1) from Definition 2.1 holds along with:

  1. the functions satisfy


    with φ(0) = φ0 a.e. in Ω.

  2. for a.e. t ∈ (0, T) and for all ΦH1, χ H01 and ζH1, (2.7b) and (2.7d) are satisfied along with

    0=Ω(νv(μ+χσ)φ)Φpdiv(Φ)dx, (2.20a)
    0=Ω1ν(p(μ+χσ)φ)χΓv(φ,σ)χdx. (2.20b)

    We say that (φ, μ, σ, v, p) is a weak solution to the CHD system with the logarithmic potential ψlog if property (c2) in Definition 2.1 holds along with properties (g) and (h).

Remark 2.2

The variational equality (2.20a) comes naturally from (2.7a) when we neglect the viscosities η(φ) and λ(φ). Meanwhile, the variational equality (2.20b) arises from the weak formulation of the elliptic problem obtained from taking the divergence of Darcy’s law (2.18) in conjunction with the equation (1.1a) and the boundary condition (2.19).

Theorem 3

Under (A1)-(A3), (B1) and (B2), there exists a weak solution (φ, μ, σ, v, p) to the CHD model with double obstacle potential ψdo in the sense of Definition 2.3. If, in addition, (C1)-(C3) hold, then there exists a weak solution (φ, μ, σ, v, p) to the CHD model with logarithmic potential ψlog in the sense of Definition 2.3, which satisfies, for a.e. t ∈ (0, T), the inequality (2.15) with η(φ) ≡ 0. Furthermore, in both cases, 0 ≤ σ ≤ 1 a.e. in ΩT.

3 Proof of Theorem 1 – Time dependent solutions

The standard procedure is to approximate the singular potentials with a sequence of regular potentials, employ Lemma 2.2 to obtain approximate solutions, derive uniform estimates and pass to the limit.

3.1 Approximation potentials and their properties

3.1.1 Double obstacle potential

We point out that in order to use Lemma 2.2 the approximate potential should at least belong to C2(ℝ). Fix δ > 0, which serves as the regularisation parameter, and we define

β^do,δ(r)=12δr1+δ22+δ24 for r1+δ,16δ2(r1)3 for r(1,1+δ),0 for |r|1,16δ2(r+1)3 for r(1δ,1),12δr+1+δ22+δ24 for r1δ. (3.1)

Formally, it is easy to see that β̂do,δ(r) → 𝕀[−1,1](r) pointwise as δ → 0, and so

ψdo,δ(r):=β^do,δ(r)+12(1r2) (3.2)

will serve as our approximation for the double obstacle potential. Let βdo,δ(r) = β̂do,δ(r) = (r + ψdo,δ(r)) ∈ C1(ℝ) denote the derivative of the convex part β̂do,δ, which has the form

βdo,δ(r)=1δr1+δ2 for r1+δ,12δ2(r1)2 for r(1,1+δ),0 for |r|1,12δ2(r+1)2 for r(1δ,1),1δr+1+δ2 for r1δ. (3.3)

Then, it is clear that βdo,δ is Lipschitz continuous with 0βdo,δ(r)1δ for all r ∈ ℝ. We now state several useful properties regarding the approximations β̂do,δ and ψdo,δ, whilst deferring the proof to the appendix.

Proposition 3.1

Let β̂do,δ and ψdo,δ be defined as above. Then, there exist positive constants C0 and C1 such that for all r ∈ ℝ and for all δ ∈ (0, 1/4),

ψdo,δ(r)C0|r|2C1, (3.4a)
δ(βdo,δ(r))22β^do,δ(r)δ(βdo,δ(r))2+1, (3.4b)
δ(βdo,δ(r))2βdo,δ(r). (3.4c)

Aside from approximating the singular potential, it would be necessary to extend the source functions bv, bφ, fv and fφ from [−1, 1] to the whole real line. Since the solution variable φ is supported in [−1, 1] (see (c1) of Definition (2.1)), the particular form of extensions outside [−1, 1] does not play a significant role and we have the flexibility to choose extensions that would easily lead to uniform estimates. Hence, unless stated otherwise we assume that bv, bφ, fv and fφ can be extended to ℝ such that fφC0(ℝ) ∩ L(ℝ), fvC0,1(ℝ) ∩ W1,∞(ℝ), bφC0(ℝ) ∩ L(ℝ), bvC0,1(ℝ) ∩ W1,∞(ℝ), and fulfil

bv(r)=0,bφ(r)=0|r|1, (3.5)
r(fφ(r)fv(r)r)<0|r|>1. (3.6)

The latter implies fφ(r) − fv(r)r is strictly negative (resp. positive) for r > 1 (resp. r < −1). For the functions stated in Example 2.1, we can consider following extensions: For r ∈ ℝ, we set



fv(r)=Aαr1,Aαr|r|1,Aαe(r+1)r1, if α<0,fv(r)=Aαe(1r)r1,Aαr|r|1,Aαr1, if α>0.

It is clear that (3.5) is fulfilled. For α = 0, we see that fv(r) = 0 and so fφ(r) − fv(r)r = fφ(r) satisfies (3.6). For α > 0, if r ≤ −1 we see that fφ(r) − fv(r)r = A(ρSα r) > 0 and if r ≥ 1 we see that fφ(r) − fv(r)r attains its maximum at r = 1 with value A(αρS) < 0. Similarly, for α < 0, if r ≥ 1 we see that fφ(r) − fv(r)r = A(α rρS) < 0 and if r ≤ −1 we see that fφ(r) − fv(r)r attains its minimum at r = −1 with value A(ρSα) > 0. Hence, the extensions fulfil (3.6). Then, we can employ Lemma 2.2 to deduce the existence of a quintuple (φδ, μδ, σδ, vδ, pδ) to (1.1)-(1.2) with ψdo,δ replacing ψ′ and source terms bv, bφ, fv and fφ modified as above. Uniform estimates will be derived in the next section and then we can pass to the limit δ → 0 to infer the existence of a weak solution to (1.1) with the double obstacle potential in the sense of Definition 2.1.

3.1.2 Logarithmic potential

We define the convex part of ψlog and its corresponding derivative as


For fixed δ > 0, the approximation of ψlog is the standard one:

ψlog,δ(r)=ψlog(1δ)+ψlog(1δ)(r(1δ))+12ψlog(1δ)(r(1δ))2 for r1δ,ψlog(r) for |r|1δ,ψlog(δ1)+ψlog(δ1)(r(δ1))+12ψlog(δ1)(r(δ1))2 for r1+δ, (3.7)

with convex part and corresponding derivative

β^log,δ(r):=ψlog,δ(r)θc2(1r2),βlog,δ(r)=β^log,δ(r)rR. (3.8)

Analogous to Proposition 3.1, we state several useful properties for β̂log,δ and ψlog,δ, whilst deferring the proof to the appendix.

Proposition 3.2

Let β̂log,δ and ψlog,δ be defined as above. Then, there exist positive constants C0, …, C3, such that for all r ∈ ℝ and for all 0 < δ ≤ min(1,θ/(4 θc)),

ψlog,δ(r)C0|r|2C1, (3.9a)
4δθ1β^log,δ(r)(|r|1)+2:=(max(0,|r|1))2, (3.9b)
δ(βlog,δ(r))22θβ^log,δ(r)+C2C3δ(βlog,δ(r))2+1, (3.9c)
δ(βlog,δ(r))2θβlog,δ(r). (3.9d)

Once again, we extend the source functions bv and bφ from [−1, 1] to the whole real line so that (3.5) is satisfied, and instead of (3.6) we consider extensions of fv and fφ from [−1, 1] to the whole real line that satisfy

rfv(r)fφ(r)>0 for r[1,r0],=0 for |r|2r0,<0 for r[r0,1], (3.10)

with continuous interpolation in [r0, 2r0] and [−2r0, −r0] for some fixed constant r0 ∈ (1, 2). For instance, with the functions fv and fφ introduced in Example 2.1, we can take as extensions

fv(r)=αA1r for |r|1,αAr for |r|1,fφ(r)=αA for r2r0,A(ρSα)2r01(r2r0)αA for r[1,2r0],ρSAr for |r|1,A(ρS+α)2r01(r+2r0)αA for r[2r0,1],αAfor r2r0,

with r0=32. Then, using ρSα > 0 it is clear that r fv(r) − fφ(r) fulfils (3.10).

3.2 Existence of approximate solutions

To unify our analysis, we use the notation

ψδ=ψdo,δ for ψdo,ψlog,δ for ψlog,Θc=1 for ψdo,θc for ψlog,

and denote by β̂δ the convex part of ψδ. Furthermore, we emphasise that the property (3.6) is assumed when we work with the double obstacle potential ψdo, and (3.10) is assumed when we work with the logarithmic potential ψlog. Due to Proposition 3.1 and 3.2 and by the definition of ψδ(⋅), we can employ Lemma 2.2, and infer, for every δ ∈ (0, 1), the existence of a weak solution quintuple (φδ, μδ, σδ, vδ, pδ) to (1.1)-(1.2) with ψδ in (1.1d). More precisely,

div(vδ)=Γv(φδ,σδ)a.e. in ΩT, (3.11a)
μδ=ψδ(φδ)Δφδχσδa.e. in ΩT, (3.11b)


0=ΩT(vδ,pδ):Φ+νvδΦ(μδ+χσδ)φδΦdx, (3.12a)
0=tφδ,ζH1+Ωμδζ+(φδvδ+φδΓv(φδ,σδ)Γφ(φδ,σδ))ζdx, (3.12b)
0=Ωσδζ+h(φδ)σδζdx+ΣK(σδ1)ζdH, (3.12c)

for a.e. t ∈ (0, T) and for all ΦH1 and ζH1.

3.3 Uniform estimates

We begin with an auxiliary result for the product of βlog,δ and the source terms, whose proof can be found in the appendix.

Lemma 3.1

Let βlog,δ denote the derivative of (3.8), and let bvC0,1(ℝ) ∩ W1,∞(ℝ), bφC0(ℝ) ∩ L(ℝ), fvC0,1(ℝ) ∩ W1,∞(ℝ), and fφC0(ℝ) ∩ L(ℝ) be given such that (C1), (C2), (3.5) and (3.10) are satisfied. Then, there exists δ0 > 0 and a positive constant C independent of δ ∈ (0, δ0) such that for all δ ∈ (0, δ0), s ∈ ℝ and r ∈ ℝ,

(rΓv(r,s)Γφ(r,s))βlog,δ(r)C(1+|s|+|r|). (3.13)

In the following, we derive uniform estimates in δ, and denote by C a generic constant independent of δ which may change its value even within one line.

3.3.1 Nutrient estimates

Choosing ζ = σδ in (3.12c) and using the non-negativity of h(⋅) leads to

Ω|σδ|2dx+KΣ|σδ|2dHKΩσδdHK2σδL2(Σ)2+K2|Σ|, (3.14)

from which we deduce that σδ is uniformly bounded in L(0, T; H1). Elliptic regularity additionally yields

σδL(0,T;H2)C. (3.15)

By a standard argument with the comparison principle, testing with ζ = − (σδ) and ζ = (σδ - 1)+ will yield the pointwise a.e. boundedness of σδ:

0σδ1a.e. in ΩT.

Then, the continuous embedding H2L and the assumptions on bv, bφ, fv and fφ lead to

ΓφL(0,T;L)+ΓvL(0,T;L)C. (3.16)

3.3.2 Estimates from energy identity

Thanks to (3.16), the function u : = 𝓓(Γv) ∈ H1 satisfies the estimate

uW1,pCΓvLpCp(1,). (3.17)

Technically, we should stress the dependence of u on δ, but in light of the uniform estimate (3.17) we infer that uδ is bounded in L(0, T; W1,p) for any p ∈ (1, ∞). Henceforth, we drop the index δ and reuse the variable u.

Choosing Φ = vδu in (3.12a), ζ = μδ + χσδ in (3.12b), using (3.11b) and summing the resulting identities, we obtain

ddtΩψδ(φδ)+12|φδ|2dx+Ω|μδ|2+2η(φδ)|Dvδ|2+ν|vδ|2dx=Ωχμδσδ+(ΓφφδΓv)(βδ(φδ)ΘcφδΔφδ)dx+Ω2η(φδ)Dvδ:Du+νvδudxΩ(βδ(φδ)ΘcφδΔφδ)φδudx, (3.18)

where we used the fact that ψδ is a quadratic perturbation of a convex function β̂δ, in conjunction with [46, Lem. 4.1] to obtain the time derivative of the energy. By Young’s inequality and the estimates (3.15), (3.16) and (3.17), we find that

Ω2η(φδ)Dvδ:Du+νvδuχμδσδdx+Ω(ΓφφδΓvφδu)(Θcφδ+Δφδ)dxη(φδ)DvδL22+ν2vδL22+14μδL22+2εΔφδL22+C1+φδH12 (3.19)

for a positive constant ε yet to be determined. It remains to control the two terms with βδ(φδ). Integrating by parts and employing the definition of u = 𝓓(Γv) leads to

Ωβδ(φδ)φδudx=Ω(β^δ(φδ))udx=1|Σ|ΩΓvdxΣβ^δ(φδ)dHΩβ^δ(φδ)Γvdx. (3.20)

Using (3.17) and the relation β^δ(r)=ψδ(r)Θc2(1r2), we obtain

Ωβ^δ(φδ)Γvdx=Ωψδ(φδ)Θc2(1φδ2)ΓvdxC1+ψδ(φδ)L1+φδL22. (3.21)

Meanwhile, using (3.4b), (3.4c), (3.9c), (3.9d) and the trace theorem, it follows that

β^δ(φδ)L1(Σ)C+Cδβδ(φδ)L2(Σ)2Cδβδ(φδ)L22+δβδ(φδ)φδL22+CCβ^δ(φδ)L1+Ωβδ(φδ)|φδ|2dx+CC1+φδL22+ψδ(φδ)L1+CΩβδ(φδ)|φδ|2dx (3.22)

for positive constants C and C* independent of δ. To deal with the remaining term we use the definition (3.3) of βdo,δ and (3.5)-(3.6) to deduce that, for any s ≥ 0 and r ∈ ℝ,

βdo,δ(r)(Γφ(r,s)rΓv(r,s))=0 for |r|1 and s0,<0 for |r|>1 and s0,

which implies that

Ω(Γφ(φδ,σδ)φδΓv(φδ,σδ))βdo,δ(φδ)dx0. (3.23)

Meanwhile, for the other part with βlog,δ, we use (3.13) and (3.15) to obtain

Ω(Γφ(φδ,σδ)φδΓv(φδ,σδ))βlog,δ(φδ)dxC(1+φδL1). (3.24)

Combining (3.23) with (3.24), we find that


and so, when substituting (3.19)-(3.22) and (3.24) into (3.18), we arrive at


Testing (3.11b) with − AΔφδ for some positive constant A, integrating by parts and using n φδ = 0 on Σ and (3.15) yields

AΩ|Δφδ|2+βδ(φδ)|φδ|2dx=AΩ(μδ+χσδ)φδ+Θc|φδ|2dxC1+φδL22+14μδL22. (3.25)

Then, summing up the last two inequalities and choosing A > C* and ε < A4 yields

ddtψδ(φδ)L1+φδL22Cψδ(φδ)L1+φδL22+μδL22+(βδ(φδ))1/2φδL22+ΔφδL22+(η(φδ))1/2Dvδ2+νvδL22C. (3.26)

Before applying a Gronwall argument, we note that for the double obstacle potential, the assumption (B2) implies β̂do,δ(φ0) = 0, and for the logarithmic potential, using the property β̂log,δ(s) ≤ β̂log(s) for s ∈ [−1, 1] together with (C3) we infer β̂log,δ(φ0) is uniformly bounded. Hence, for 0 < δ < min (1, θ/(4 θc), δ0) = : δ*, we see that


Integrating (3.26) in time from 0 to s ∈ (0, T], using (3.4a), (3.9a), Korn’s inequality and elliptic regularity theory, we deduce the uniform estimate

ψδ(φδ)L(0,T;L1)+φδL(0,T;H1)L4(0,T;H2)+μδL2(0,T;L2)+(βδ(φδ))1/2φδL2(0,T;L2)+vδL2(0,T;H1)C, (3.27)

where the L4(0, T; H2)-regularity comes from revisiting the equality in (3.25) and employing the boundedness of ∇ μδ in L2(QT) and of φδ in L(0, T; H1(Ω)), the non-negativity of βδ , leading to


and also elliptic regularity theory. By the Sobolev embedding H1L6 and (3.27), it follows that


A similar argument together with (3.16) shows that φδ Γv is bounded in L2(0, T; L32 ). Then, from (2.7b) we obtain

tφδL2(0,T;(H1))+div(φδvδ)L2(0,T;L3/2)C. (3.28)

Furthermore, we find that the mean value (φδ)Ω satisfies


and so

(φδ)ΩH1(0,T)C. (3.29)

In particular, by the fundamental theorem of calculus, it holds that

(φδ)Ω(r)(φδ)Ω(s)=srt(φδ)Ω(t)dtC|rs|12. (3.30)

3.4 Estimates for the mean value of the chemical potential

In order to pass to the limit δ → 0 rigorously, it remains to derive uniform estimates for μδ, βδ(φδ) and pδ in L2(0, T; L2). To do so we appeal to the method introduced in [33], which involves first deducing that the limit φ of φδ has mean value in the open interval (−1, 1) for all times. We first state a useful auxiliary result, whose proof is deferred to the appendix.

Proposition 3.3

For δ ∈ (0, 1), let βlog,δ denote the derivative of (3.8). Then, there exist positive constants c1 and c2 independent of δ such that

rβlog,δ(r)|βlog,δ(r)|c1|r|c2rR. (3.31)

Remark 3.1

This is a different type of estimate to the more commonly used inequality (see e.g. [13, (2.12)]):

rβlog,δ(r)c~0|βlog,δ(r)|c~1rR (3.32)

with positive constant 0 and non-negative constant 1 that are independent of δ, provided δ is sufficiently small. We found that (3.32) is insufficient for our analysis as the constant 0 is not quantified.

Now, using reflexive weak compactness arguments (Aubin–Lions theorem) and [50, Sec. 8, Cor. 4], for δ → 0 along a non-relabelled subsequence, we infer that

φδφweakly- in H1(0,T;(H1))L(0,T;H1)L2(0,T;H2),φδφstrongly in C0([0,T];Lr)L2(0,T;W1,r) and a.e. in ΩT,σδσweakly- in L(0,T;H2),μδξweakly in L2(0,T;L2),vδvweakly in L2(0,T;H1),div(φδvδ)θweakly in L2(0,T;L32), (3.33)

for some limit functions ξL2(0, T; L2), θL2(0, T; L32 ) and for all r ∈ [1, 6). The interpolation inequality fH1CfL21/2fH21/2, the boundedness of φδ - φ in L2(0, T; H2) and the strong convergence φδφ in L(0, T; L2) also allow us to deduce that

φδφ strongly in L4(0,T;H1). (3.34)

Let λL4(0, T; L3) be an arbitrary test function, then (3.34) implies λ φδλ φ strongly in L2(0, T; L2) and λφδλφ strongly in L2(0, T; L65 ). Using the weak convergence of vδv in L2(0, T; H1) and the product of weak-strong convergence we obtain

0TΩdiv(φδvδ)λdxdt0TΩdiv(φv)λdxdtas δ0. (3.35)

This implies div(φδvδ) → div(φv) weakly in L43(0,T;L32) as δ → 0. Since L2(0,T;L32)L43(0,T;L32), by uniqueness of weak limits we obtain div(φv) = θ. Using the assumption on Γφ and the above convergences, we can pass to the limit in (2.7b) to obtain

tφ,ζH1+Ωdiv(φv)ζdx=Ωξζ+Γφ(φ,σ)ζdx (3.36)

for a.e. t ∈ (0, T) and for all ζH1. Technically, one would multiply (2.7b) with a function κ Cc (0, T), integrate the resulting product over (0, T), pass to the limit δ → 0 and then recover (3.36) with the fundamental lemma of calculus of variations.

Now, for the obstacle potential, the uniform boundedness of ψδ(φδ) in L1(0, T; L1) and (3.4b) imply δβdo,δ(φδ) is uniformly bounded in L2(0, T; L2), and so δ βdo,δ(φδ) → 0 strongly in L2(0, T; L2). However, from the definition of βdo,δ we have the pointwise convergence

δβdo,δ(r)g(r):=r1 if r1,0 if |r|1,r+1 if r1,as δ0, (3.37)

which implies that (c.f. [8, Proof of Thm. 2.2])

|φ|1a.e. in ΩT. (3.38)

For the logarithmic potential, we use (3.9b) and the uniform boundedness of β̂log,δ(φδ) in L1(0, T; L1) to obtain that


Since φδφ a.e. in ΩT and strongly in L2(0, T; L2), passing to the limit δ → 0 in the inequality above also implies (3.38). From this we claim that φΩ(t) ∈ (−1, 1) for all t ∈ [0, T]. Indeed, choosing ζ = 1 in (3.36) leads to

tφ,1H1+Ωφvdx=ΩΓφ(φ,σ)φΓv(φ,σ)dxfor a.e.t(0,T). (3.39)

Suppose to the contrary there exists a time t* ∈ (0, T] such that φΩ(t*) = 1 and (3.39) holds. Due to (3.38), this implies φ(t*, x) ≡ 1 a.e. in Ω and thus ∇ φ(t*, x) ≡ 0 a.e. in Ω. Using (3.39) and (3.5)-(3.6), we obtain


Hence, by continuity of t ↦ (φ(t))Ω, the mean value (φ(t))Ω must be strictly decreasing in a neighbourhood of t*, i.e., (φ(t))Ω > 1 for t < t*, which contradicts (3.38). Using a similar argument and the assumption fφ(−1) + fv(−1) > 0 leads to the conclusion that (φ(t))Ω > − 1 for all t ∈ (0, T]. In particular, together with (B2), we have (φ(t))Ω ∈ (−1, 1) for all t ∈ [0, T].

This allows us to derive uniform estimates on the mean value of μδ. Integrating (3.11b) and taking the modulus on both sides gives

Ωμδ(t)dxΩ|βδ(φδ(t))|+Θc|φδ(t)|+χ|σδ(t)|dx (3.40)

for a.e. t ∈ (0, T). Using (3.31) and the fact

|βdo,δ(r)|rβdo,δ(r) for all rR (3.41)

(which unfortunately does not hold for βlog,δ, hence the necessity of Proposition 3.3), we deduce that

|βδ(r)|rβδ(r)+c1|r|+c2 for all rR,

and so, we arrive at

Ωμδ(t)dxΩ|βδ(φδ(t))|+Θc|φδ(t)|+χ|σδ(t)|dx,Ωφδ(t)βδ(φδ(t))+(c1+Θc)|φδ(t)|+χ|σδ(t)|+c2dx. (3.42)

Together with the following equality obtained from testing (3.11b) with ζ = φδ,


we see that

Ωμδ(t)dxΩμδ(t)φδ(t)dx+C1+φδ(t)L22+σδ(t)L22. (3.43)

Now, let fδHn2L02 be the unique solution to the Neumann problem

Δfδ=φδ(t)(φδ(t))Ω in Ω,nfδ=0 on Σ, (3.44)

where by Poincaré’s inequality, we have

fδH1Cφδ(t)L2. (3.45)

Testing (3.12b) with fδ, integrating by parts and rearranging yields


Plugging in this identity into (3.43) and rearranging again, we deduce that

1|(φ(t))Ω|supt(0,T)|(φδ(t)φ(t))Ω|Ωμδ(t)dxCσδ(t)L22+φδ(t)L22tφδ(t),fδH1Ω(div(φδ(t)vδ(t))Γφ(φδ(t),σδ(t)))fδdx (3.46)

for a.e. t ∈ (0, T). Recalling (3.29)-(3.30), we have the equiboundedness and equicontinuity of {(φδ)Ω}δ∈(0,1) so that by the Arzelá–Ascoli theorem,

(φδ(t))Ω(φ(t))Ωstrongly in C0([0,T])as δ0.

along a non-relabelled subsequence. Then, one can find an index δ1 ∈ (0, 1) such that for all δ < min(δ1, δ*) = : δ2, where δ* is defined after (3.26),


where we know that the minimum is positive thanks to the property |(φ(t))Ω| < 1 for all t ∈ [0, T]. Then, we see that the prefactor on the left-hand side of (3.46) is bounded below by 12 min[0,T] (1 − |(φ(t))Ω|), which is positive uniformly in t. As the right-hand side of (3.46) is uniformly bounded in L2(0, T) by previously established uniform estimates, we obtain {(μδ)Ω}δ∈(0,δ2) is bounded in L2(0, T), and the Poincaré inequality gives

μδL2(0,T;L2)C. (3.47)

Let us mention that if instead of (3.31) we employ the estimate (3.32), then we arrive at


and ultimately


Since 0 is not quantified, we may not be able to rule out the situation where 0 < 1, which may imply that the prefactor 1|(φ(t))Ω|min(1,c~0) is negative.

The uniform estimate (3.47) for μδ allows us to infer further estimates for βδ(φδ) and pδ. Indeed, testing (3.11b) with βδ(φδ) yields


Integrating this identity in time from 0 to T, using the non-negativity of βδ (⋅), (3.27) and (3.47), it follows that

βδ(φδ)L2(0,T;L2)C. (3.48)

For the pressure pδ, we invoke Lemma 2.1 to deduce the existence of qδ := 𝓓(pδ) ∈ H1 such that for a positive constant C depending only on Ω,

qδH1CpδL2. (3.49)

Then, testing (3.12a) with Φ = qδ yields


Applying Young’s inequality and using the uniform estimates (3.15), (3.27), (3.47) and (3.49) leads to

pδL2(0,T;L2)C. (3.50)

Following similar arguments in [16], with the uniform boundedness of μδ + χ σδ + Θc φδ in L2(0, T; Lp) for p = 6 if d = 3 and for any p ∈ [2, ∞) if d = 2, we can deduce from (3.11b) further uniform estimates


3.5 Passing to the limit

Let us first consider the double obstacle case. In addition to the convergence statements in (3.33), we further obtain

μδμweaklyin L2(0,T;H1),βdo,δ(φδ)τweaklyin L2(0,T;L2),pδpweaklyin L2(0,T;L2),

for some limit function τL2(0, T; L2). Moreover, due to (3.47) we have ξ = ∇μ, which allows us to fully recover (2.7b) in the limit. To obtain (2.7a) and (2.7d) in the limit we refer the reader to similar arguments as outlined in [19], while the divergence equation (1.1a) can be obtained by similar arguments in [20]. It remains to show (2.11) is recovered in the limit δ → 0 from (3.11b). By arguing as in [31, Sec. 5.2], using the weak convergence βdo,δ(φδ) ⇀ τ in L2(0, T; L2), the strong convergence φδφ in L2(0, T; L2), and the maximal monotonicity of the subdifferential 𝕀[−1,1] we can infer that τ is an element of the set 𝕀[−1,1](φ), which implies that for any ζ ∈ 𝓚 and a.e. t ∈ (0, T),


Hence, testing (3.11b) (where ψδ = ψdo,δ) with ζφ and passing to the limit δ → 0 allows us to recover (2.11). This completes the proof of Theorem 1 for the double obstacle potential.

For the logarithmic case, the uniform estimate for βlog,δ(φδ) in L2(0, T; L2) allow us to infer, by the arguments in [13, Sec. 4] or [35, Sec. 3.3], that the limit φ satisfies the tighter bounds

|φ(x,t)|<1a.e. in ΩT.

Furthermore, by the a.e. convergence of φδ to φ we have βlog,δ(φδ) → βlog(φ) a.e. in ΩT. Hence, passing to the limit δ → 0 in (3.11b) we can recover (2.12). Meanwhile, the inequality (2.15) is obtained from integrating (3.18) over (0, t) for t ∈ (0, T) and then passing to the limit with the compactness assertions (3.33), weak lower semicontinuity, and Fatou’s lemma. This completes the proof of Theorem 1 for the logarithmic potential.

4 Proof of Theorem 2 – Stationary solutions

As with the time-dependent case, we extend bv, bφ, fv and fφ from [− 1, 1] to ℝ such that fφC0(ℝ) ∩ L(ℝ), fvC0,1(ℝ) ∩ W1,∞(ℝ), bφC0(ℝ) ∩ L(ℝ), bvC0,1(ℝ) ∩ W1,∞(ℝ), and fulfil (3.5), (3.6) (if ψ = ψdo) or (3.10) (if ψ = ψlog).

4.1 Approximation scheme

We consider a smooth function ĝ :ℝ → [0, 1] such that ĝ(r) = 1 for r ≥ 3 and ĝ(r) = 0 for r ≤ 2, and define F : L2(Ω) → ℝ as

F(v):=CFg^(1|Ω|vL22) for vL2(Ω),

where CF is a positive constant to be specified later. Furthermore, we denote by y(r, s) the function


and introduce, for δ ∈ (0, 1), the cut-off operator 𝓓δC0,1(ℝ) defined as

Tδ(s)=1δ if s1δ,s if |s|1δ,1+δ if s1+δ. (4.1)

It is clear that 𝓓δ and its derivative Tδ are bounded independently of δ. Then, we seek a solution φδW := Hn2 to the approximate system

δβδ(φδ)+F(φδ)φδ+Δ(Δφδψδ(φδ)+χσδ)=y(φδ,σδ)vδTδ(φδ) in Ω,nφδ=n(Δφδ+χσφδ)=0 on Σ, (4.2)

where σδH2 is the unique non-negative and bounded solution to the nutrient subsystem

0=Δσδh(φδ)σδ in Ω,nσδ=K(1σδ) on Σ, (4.3)

and vδ is the first component of the unique solution (vδ, pδ) to the Brinkman subsystem

div(T(vδ,pδ))+νvδ=(ψδ(φδ)Δφδ)Tδ(φδ) in Ω,div(vδ)=Γv(φδ,σδ) in Ω,T(vδ,pδ)n=(2η(φδ)Dvδ+λ(φδ)div(vδ)IpδI)n=0 on Σ. (4.4)

We aim to use pseudomonotone operator theory, akin to the methodology used in [33], to deduce the existence of at least one solution φδW to (4.2) for each δ ∈ (0, 1). Then, we derive enough uniform estimates to pass to the limit δ → 0 in order to prove Theorem 2. The new element in the analysis is how we treat the convection term.

4.2 Preparatory results

For fixed uW := Hn2 it is clear that there is a unique solution σuH2 to the nutrient subsystem (4.3), and

σu1σu2H2Cu1u2L2 (4.5)

for any pair u1, u2W with corresponding solutions σu1, σu2. For the Brinkman subsystem (4.4), invoking Lemma 2.3 with c = uW, g = Γv(u, σu) ∈ H1, f = (ψδ(u) − Δ u) ∇ 𝓓δ(u) ∈ L54 , and h = 0, we can immediately deduce for any uW and δ ∈ (0, 1) the existence of a unique strong solution (vu,pu)W2,54×W1,54.

Although this is sufficient for the upcoming analysis, we mention that an analogous continuous dependence result to (4.5) for the Brinkman subsystem (4.4) is not available from Lemma 2.3, and thus for completeness we provide the following result under a modification of the cut-off operator 𝓓δ.

Proposition 4.1

Suppose 𝓓δC1,1(ℝ), and let (v1, p1), (v2, p2) ∈ W2,54×W1,54 be the solutions to (4.4) corresponding to data u1, u2W, respectively. Then, there exists a positive constant C independent of the differences û = u1u2, = v1v2 and = p1p2 such that


One example of the modified cut-off operator 𝓓δC1,1(ℝ) is the function

132δ if s1δ,132δ12δ(s(1δ))2 if 12δs1δ,s if |s|12δ,32δ1+12δ((δ1)s)2 if 1+δs1+2δ,1+32δ if s1+δ.


For convenience we write η^=η(u1)η(u2),ψ^δ=ψδ(u1)ψδ(u2),Tδ^=Tδ(u1)Tδ(u2) and Γ̂v := Γv(u1, σu1) − Γv(u2, σu2). Let ŵ := 𝓓(Γ̂v) ∈ H1 where by (4.5) it holds that

wH1CΓ^vL2Cu^L2. (4.6)

Then, testing the difference of (4.4)1 with ŵ yields


The left-hand side can be bounded below by

LHSν2v^L22+η02Dv^L22ν2w^L22CDw^L22C2Dv2L22u^L2, (4.7)

while the right-hand side can be bounded above by


In light of the regularities u1, u2W and v1, v2 W2,54 , the polynomial growth of ψδ leading to the following difference estimate


the Lipschitz continuity and boundedness of Tδ , as well as (4.6), we find that


for some small constant ε > 0. Invoking Korn’s inequality and (4.6) in (4.7) and combining with the above estimate gives


For the pressure, we test the difference of (4.4)1 with := 𝓓() ∈ H1 which satisfies


and in turn we obtain


This completes the proof. □

The unique solvability of the nutrient and Brinkman subsystems in turn yields the following result.

Lemma 4.1

For each δ ∈ (0, 1), the operator 𝓐 : WW* defined as


is strongly continuous, i.e., if unu in W, then 𝓐un → 𝓐u in W*.


Let {un}n∈ℕW be a sequence of functions such that unu in W. Denoting by σn and (vn, pn) the corresponding unique solutions to the nutrient subsystem (4.3) and Brinkman subsystem (4.4) where φδ = un. Then, we easily infer that

σnH2C,0σn1 a.e. in Ω,

for a positive constant C independent of n. Furthermore, using the assumptions on Γv by Lemma 2.1, there exists a function un := 𝓓(Γv(un, σn)) ∈ H1 satisfying


with C independent of δ. Testing (4.4)1 (for φδ = un) with vnun, using (4.4)2-(4.4)3 and Korn’s inequality, we obtain


Similarly, let qn := 𝓓(pn) ∈ H1 which satisfies


Testing (4.4)1 (for φδ = un) with qn and using (4.4)2-(4.4)3, we obtain


In particular, this shows that

vnH1+pnL2C1+unL2ψδ(un)L3+ΔunL2unL3Cδ, (4.8)

for a positive constant Cδ depending on δ but independent of n, thanks to the boundedness of {un}n∈ℕ in W. Hence, for fixed δ ∈ (0, 1), there exist functions σuH2, vuH1 and puL2, such that along a non-relabelled subsequence,

σnσu in H2,vnvu in H1,pnpu in L2.

It is clear that σu is the unique solution to (4.3) corresponding to u, while (vu, pu) is the unique weak solution to (4.4) corresponding to u. Indeed, for the highest order term in (4.4), we employ the dominated convergence theorem to deduce that Tδ (un) → Tδ (u) strongly in L152 . Together with Δ unΔ u in L2 and ∇ un → ∇ u in L5, we obtain

ΔunTδ(un)ΔuTδ(u) in L6/5,

while the other terms in the Brinkman system can be handled using similar compactness assertions.

By Rellich’s theorem, it is easy to see that


for all ζW. Meanwhile, integrating by parts and using div (vn) = Γv(un, σn) = : Γvn yields


on account of 𝓓δ(un) → 𝓓δ(u) in L2 and L2(Σ), vnvu in L3 and L2(Σ) (see, e.g., [19, Lem. 1.3]), and ΓvnΓvu in L2 . This shows that 𝓐 is strongly continuous. □

4.3 Existence of approximate solutions

In this section we fix δ ∈ (0, 1), and define operators A1, A2 : WW* by


Then, φδW is a weak solution to (4.2) if 〈(A1 + A2)φδ, ζW = 0 for all ζW.

Since βδ is bounded and βδ has sublinear growth, we deduce that the operator A1 is monotone and hemicontinuous. On the other hand, Lemma 4.1 together with the continuity and sublinear growth of ψδ , and the continuity and boundedness of F imply that A2 is strongly continuous. Then, by [52, Thm. 27.6] the sum A = A1 + A2 is a pseudomonotone operator.

We now claim that A is additionally coercive over W, i.e.,


Let us first treat the velocity term in 〈Au, uW. By the definition of 𝓓δ in (4.1), we see that u Tδ (u) ∇ u = 𝓓δ(u) Tδ (u) ∇ u = 𝓓δ(u) ∇ 𝓓δ(u) a.e. in Ω. Hence, by the trace theorem we have

Ωu(vuTδ(u))dx=Ω12vu|Tδ(u)|2dx=Σ12|Tδ(u)|2vundHΩ12Γv(u,σu)|Tδ(u)|2dx12vuL1(Σ)+CC(1+vuH1). (4.9)

Let uu := 𝓓(Γvu(u, σu)) ∈ H1 which satisfies


where the boundedness of Γv(u, σu) can be deduced using similar arguments in Section 3.3. Testing (4.4) for (vu, pu) with vuuu then yields


By Korn’s inequality, Hölder’s inequality and Young’s inequality we obtain


where C* is the constant in (4.9). Let us note that as a consequence of elliptic regularity and integration by parts, for uW we obtain the following useful inequalities:

uL2ΔuL21/2uL21/2,uH2CΩ(ΔuL2+uL2). (4.10)

Employing the Gagliardo–Nirenberg inequality and the elliptic estimate (4.10), we have


Meanwhile, we observe that

Tδ(s)βδ(s)=Tδ(s)βlog(s) for the logarithmic potential,0 for the obstacle potential,

since βdo,δ(s) = 0, βlog,δ(s) = βlog(s) for s ∈ [− 1 + δ, 1 − δ]. Hence, as 0 ≤ Tδ (s) ≤ 1, we see that

Ω|ψδ(u)Tδ(u)u|54dx={|u|1δ}|(βlog(u)Θcu)u|54dx{|u|1δ}|βlog(u)|54|βlog(u)|58|βlog(u)|58|u|54dx+C(1+uL22)CΩ|βlog,δ(u)|58|u|54dx+C(1+uL22)12Ωβlog,δ(u)|u|2dx+18ΔuL22+C(1+uL22), (4.11)

where for the second last inequality we used L’Hopital’s rule to confirm that βlog(s)(βlog(s))1/2 is continuous in [− 1, 1], hence bounded. Returning to (4.9), we deduce that

Ωu(vuTδ(u))dx14ΔuL22+12Ωβδ(u)|u|2+C1(1+uL22) (4.12)

for a positive constant C1 independent of u and δ. Now, computing 〈Au, uW gives

Au,uW=ΩF(u)|u|2+δβδ(u)u+|Δu|2+βδ(u)|u|2Θc|u|2dx+Ωuy(u,σu)+u(vuTδ(u))+χuΔσudxΩF(u)|u|2+δβδ(u)u+12|Δu|2+12βδ(u)|u|2dxC(1+uL22), (4.13)

for a positive constant C independent of u and δ. Recalling that F(u) = CF for uL22>3|Ω|, and so, in choosing CF = 2C, we arrive at


for uL223|Ω|, which in turn implies coercivity of A.

Invoking [52, Thm. 27.A], we deduce for every δ ∈ (0, 1) the existence of a weak solution φδW to (4.2). Setting

μδ=Δφδ+ψδ(φδ)χσφδ (4.14)

we see that the equation A φδ = 0 in W* implies

ΩμδΔζdx=ΩfδζdxζW, (4.15)

with right-hand side


Thanks to the regularity φδW, vδH1, σδH2, and the sublinear growth of βδ, we easily infer that fδL2. On the other hand, choosing ζ = 1 in (4.15) implies that f L02 . Then, by arguing as in [33, Sec. 3.1], we obtain that μδW for all δ ∈ (0, 1).

4.4 Uniform estimates

From (4.2) and (4.15), the pair (φδ, μδ) ∈ W × W satisfies

0=Ω(F(φδ)φδ+δβδ(φδ)+y(φδ,σδ)+vδTδ(φδ)Δμδ)ζdx, (4.16a)
0=Ω(βδ(φδ)ΘcφδμδχσδΔφδ)ζdx (4.16b)

for all ζL2. Returning to the proof of the coercivity of the operator A, replacing u with φδ in (4.13) gives

ΩF(φδ)|φδ|2+δβδ(φδ)φδ+12βδ(φδ)|φδ|2+12|Δφδ|212CFφδL22+C. (4.17)

If φδL223|Ω|, then as before we absorb the 12CFφδL22 on the right-hand side by the F(φδ)|φδ|2 term on the left-hand side and employ the elliptic estimate (4.10) to obtain

Ωδβδ(φδ)φδ+βδ(φδ)|φδ|2dx+φδH22C. (4.18)

If φδL22 < 3 |Ω|, then adding φδL22 to both sides of (4.17) and neglecting the non-negative term F(φδ) |φδ|2 on the left-hand side yields the uniform estimate (4.18). Hence, {φδ}δ∈(0,1) is bounded in W, and along a non-relabelled subsequence, it holds that

φδφ,σδσφ in H2,

where σφ is the unique solution to the nutrient subsystem (4.3) with data φ.

Convexity of β̂δ and β̂δ(0) = 0 imply the inequality

β^δ(s)βδ(s)s for all sR,

and together with (3.4b), (3.9b) and (4.18) we deduce that


Using (3.37) for the double obstacle potential, we deduce that for both cases the limit φ satisfies

|φ|1 a.e. in Ω. (4.19)

In particular, we have

φL22|Ω|,Tδ(φδ)φ in Lp for any p<. (4.20)

The latter can be deduced from the a.e. convergence 𝓓δ(φ) → φ, the Lipschitz continuity of 𝓓δ and hence a.e. convergence of 𝓓δ(φδ) − 𝓓δ(φ) → 0, and the dominating convergence theorem. Using the norm convergence φδL22φL22, we then infer the existence of δ3 > 0 such that φδL22 ≤ 2 |Ω| for δ ∈ (0, δ3). Subsequently, F(φδ) = 0 for δ ∈ (0, δ3), and in the sequel we will neglect the term F(φδ)φδ.

Choosing ζ = − Δ μδ in (4.16b), and choosing ζ = βδ(φδ) and also ζ = - Δ φδ in (4.16a) yields after summation and integrating by parts

μδL22+δβδ(φδ)L22Ωy(φδ,σδ)(Δφδβδ(φδ))(Θcφδ+χσδ)μδdx+ΩvδTδ(φδ)(Δφδβδ(φδ))dxC+12μδL22+ΩvδTδ(φδ)(Δφδβδ(φδ))dx (4.21)

on account of the boundedness of φδ and σδ in H2, and the estimates (3.23) and (3.24) for the term y(φδ, σδ) βδ(φδ). Let uδ := 𝓓(Γvδ) ∈ H1 which satisfies


Then, testing vδuδ with the Brinkman subsystem (4.4) yields

ΩvδTδ(φδ)(βδ(φδ)Δφδ)dx2η1/2(φδ)DvδL22νvδL22=ΩΘcφδTδ(φδ)vδ+(ψδ(φδ)Δφδ)Tδ(φδ)uδdxΩ2η(φδ)Dvδ:Duδ+νvδuδdx. (4.22)

Recalling (4.9) and (4.11), when substituting (4.22) to (4.21), we have


where ε > 0 is a constant yet to be determined. Invoking Korn’s inequality and choosing ε sufficiently small, we arrive at

μδL22+δβδ(φδ)L22+vδH12C. (4.23)

In particular, we deduce that vδv in H1 to some limit velocity v. Moreover,

δβδ(φδ)L22Cδ0 as δ0,

and so δβδ(φδ)0 in L2. Then, choosing ζ = 1 in (4.16a), using that μδW and passing to the limit δ → 0 yields

0=Ωy(φ,σ)+vφdx, (4.24)

where the convergence of the convection term follows from an analogous integration by parts as in the proof of Lemma 4.1, and the convergence (4.20). Then, we infer from (2.13) and (4.24) that the limit φ has mean value φΩ ∈ (−1, 1). Indeed, substituting φ = 1 or − 1 in (4.24) leads to a contradiction on account of (2.13), and as |φ| ≤ 1 a.e. in Ω, we have that φΩ ∈ (−1, 1).

Arguing as in the time-dependent case we can derive a uniform estimate on the mean value of μδ, and consequently


where the boundedness of βlog,δ(φδ) in L2 implies the tighter bounds

|φ|<1 a.e. in Ω

in the case of the logarithmic potential.

4.5 Passing to the limit

In (4.16a) we take ζH1 and apply integration by parts to get


Passing to the limit δ → 0 then yields (2.17). Meanwhile, (2.11) or (2.12) can be recovered in the limit δ → 0 from (4.16b) in a fashion similar to the time-dependence case, as with the recovery of (1.1a) and (2.7d).

It remains to show that the weak limit p of (pδ)δ∈(0,δ3) together with v constitutes a solution to the corresponding Brinkman system in the sense of (2.7a). From the definition (4.14) of μδ we observe from (4.4) that (vδ, pδ) satisfies

0=Ω2η(φδ)Dvδ:DΦ+λ(φδ)div(vδ)div(Φ)+νvδΦdxΩpδdiv(Φ)+(μδ+χσδ)Tδ(φδ)Φdx (4.25)

for ΦH1. For the last term, after integrating by parts and using 𝓣δ(φδ) Φφ Φ in L2(Ω) and in L2(Σ), μδμ in L4(Ω) and in L2(Σ), we see that


for all ΦH1. Hence, passing to the limit δ → 0 in (4.25) allows us to recover (2.7a) and thus the quintuple (φ, μ, σ, v, p) is a stationary solution in the sense of Definition 2.2.

Moreover, from the above estimates and weak lower semicontinuity of norms, we know that

φH2+μH1+σH2+vH1+pL2C. (4.26)

Then, from (2.17) and elliptic regularity, we deduce that

μH2C. (4.27)

In light of this improved regularity and the Sobolev embedding H2L, it is easy to see that (μ + χ σ) ∇ φLq where q < ∞ for d = 2 and q = 6 for d = 3. Invoking Lemma 2.3, we then infer

vW2,q+pW1,qC, (4.28)

which completes the proof.

5 Proof of Theorem 3 – Darcy’s law

We can adapt most of the arguments and estimates from the proof of Theorem 1. The main idea is to consider a weak solution quintuple (φδ, μδ, σδ, vδ, pδ) to the CHB model (1.1)-(1.2) with stress tensor


where we have set η(⋅) = λ(⋅) = δ. Proceeding as in the proof of Theorem 1 we obtain the uniform estimates (3.15), (3.16) and

ψδ(φδ)L(0,T;L1)+φδL(0,T;H1)L2(0,T;H2)+μδL2(0,T;L2)+(βδ(φδ))1/2φδL2(0,T;L2)+vδL2(0,T;Ldiv2)+δDvδL2(0,T;L2)C, (5.1)

where in the above ψδ and βδ denote the approximations to either singular potentials and the derivatives of the corresponding convex part. From the first equality of (3.25) with A = 1, it holds that


so that by elliptic regularity we can infer

φδL4(0,T;H2)C. (5.2)

Moreover, by the Gagliardo–Nirenberg inequality, we find that


so that from (2.7b) and previous uniform estimates we arrive at

tφδL8/5(0,T;(H1))+φδvδL8/5(0,T;L6/5)C, (5.3)
tφδ+div(φδvδ)L2(0,T;(H1))+(φδ)ΩW1,8/5(0,T)C, (5.4)
|(φδ)Ω(r)(φδ)Ω(s)|C|rs|3/8r,s(0,T). (5.5)

Let us mention that the sum t φδ + div(φδ vδ) has better temporal integrability than either of its constituents, a fact which will play an important role for deriving uniform estimates for (μδ)Ω below.

By reflexive weak compactness arguments and [50, Sec. 8, Cor. 4], for δ → 0 along a non-relabelled subsequence, it holds that for any r ∈ [1, 6),

φδφ weakly-* in W1,85(0,T;(H1))L(0,T;H1)L4(0,T;H2),φδφ strongly  in C0([0,T];Lr)L4(0,T;W1,r) and a.e. in ΩT,σδσ weakly-* in L(0,T;H2),vδv weakly  in L2(0,T;L2),div(φδvδ)θ weakly  in L85(0,T;L65).

The identification θ = div(φ v) follows analogously as in Section 3.4, where the assertion (3.35) now holds for arbitrary λL4(0, T; L6) by the strong convergence ∇ φδ → ∇ φ in L4(0, T; L3) and the weak convergence vδ harpoonup v in L2(0, T; L2).

In order to obtain uniform estimates for the chemical potential μδ in L2(0, T; L2), we again follow the argument in Section 3.4. Namely, we pass to the limit δ → 0 in (2.7b) to obtain (3.36), and use the uniform boundedness of ψδ(φδ) in L1(0, T; L1) from (5.1) to obtain that the limit φ satisfies the pointwise bound (3.38). Choosing ζ = 1 in (3.36) leads to (3.39) and obtain by contradiction argument that (φ(t))Ω ∈ (−1, 1) for all t ∈ [0, T].

Defining fδ Hn2L02 as the unique solution to (3.44) which satisfies (3.45). Then, the right-hand side of (3.46) can be estimated as

RHSC(σδ(t)L22+φδ(t)L22)+Γφ(φδ(t),σδ(t))L2fδL2+tφδ(t)+div(φδ(t)vδ(t))(H1)fδH1, (5.6)

which is bounded in L2(0, T) by (5.4). This modification allows us to infer that (μδ)Ω is uniformly bounded in L2(0, T), whereas simply using (5.3) would only give the uniform boundedness of (μδ)Ω in L85 (0, T). Hence, we recover the uniform L2(0, T; L2)-estimate (3.47) for μδ and also (3.48) for βdo,δ(φδ) and βlog,δ(φδ). Moreover, setting η(⋅) = δ and η1 = λ1 = δ, we obtain as in the end of Section 3.4 the uniform L2(0, T; L2)-estimate (3.50) for pδ.

Then, proceeding as in the proof of Theorem 1, we can recover (2.7b), (2.7d) and (2.11) (resp. (2.12)) for the double obstacle (resp. logarithmic) case in the limit δ → 0, whereas recovery of (2.20a), (2.20b), the improved regularity p L85 (0, T; H1) and the boundary condition (2.19) follow from similar arguments outlined in [20, Sec. 4.2].

A Properties of approximation potentials

A. 1 Proof of Proposition 3.1


As ψdo,δ is bounded for ∣r∣ ≤ 1 + δ, it suffices to show (3.4a) holds for ∣r∣ > 1 + δ. By Young’s inequality it is clear that for r > 1 + δ with δ ∈ (0, 1/4),


and a similar assertion holds also for r < −1 − δ. This establishes (3.4a).

From the definitions of β̂do,δ and βdo,δ we see that


for r ∈ (1, 1 + δ), and


for r > 1 + δ. Similar assertions also hold for the cases r ∈ (−1 − δ, −1) and r < −1 − δ, which then yield (3.4b).

A straightforward computation shows

βdo,δ(r)=δ(βdo,δ(r))2=1δ for |r|1+δ,0 for |r|1,δ(βdo,δ(r))2=1δ3(r1)21δ2(r1)=βdo,δ(r) for r(1,1+δ),δ(βdo,δ(r))2=1δ3((r+1))21δ2(r+1)=βdo,δ(r) for r(1δ,1),

and so (3.4c) is established.□

A.2 Proof of Proposition 3.2


For r ≥ 1 − δ with δθ4θc , a short calculation shows that


where the inequality βlog(1-δ) ≥ 2 θc comes from the facts βlog(1δ)=θδ(2δ) and δ(2 − δ) ≤ 2 δ θ2θc. Then, it is easy to see that (3.9a) holds with the help of Young’s inequality. Analogously, using β̂log(δ − 1) > 0, βlog(δ − 1) = βlog(1 − δ) and βlog(δ − 1)(r - (δ − 1)) ≥ 0 for r ≤ −1 + δ, we infer that (3.9a) also holds for r ≤ −1 + δ with δθ4θc . Meanwhile for ∣r∣ ≤ 1- δ, using the non-negativity of β̂log yields


This completes the proof of (3.9a).

For ∣r∣ ≤ 1, we see that (|r|1)+2=0 and so (3.9b) holds trivially due to the non-negativity of β̂log,δ. For r > 1, the definition of β̂log,δ gives


and similar arguments yield β^log,δ(r)θ4δ(|r|1)2 for r < −1, which shows (3.9b).

For (3.9c) we first observe that δ (βlog,δ(0))2 = 0 = β̂log,δ(0) and for δ ∈ (0, 1] we have


which imply


Integrating the first inequality from 0 to r ∈ (0, 1 − δ] and the second inequality from r ∈ [−1 + δ, 0) to 0 yields


Taking note that over [−1 + δ, 1 − δ], β̂log,δ(r) is bounded uniformly in δ ∈ (0, 1], and so we easily infer the upper bound 2 θ β̂log,δ(r) ≤ C2 (δ (βlog,δ(r))2 + 1) for some positive constant C2 > 0 holding for all r ∈ [−1 + δ, 1 − δ]. Meanwhile, a direct calculation shows that for r ≥ 1 − δ and δ ∈ (0, 1] we have


on account of 1212δ1, the positivity of log((2 − δ)/δ) (r − (1 − δ)) and the boundedness of β̂log(1 − δ) and δ log ((2 − δ)/δ)2 for δ ∈ (0, 1]. An analogous calculation leads to similar inequalities for r ≤ −1 + δ, and thus (3.9c) is established.

For (3.9d) a straightforward calculation using 12δ1 gives


This completes the proof.□

A. 3 Proof of Lemma 3.1


We define


Then, due to (3.5), G(r, s) = 0 for ∣r∣ ≥ 1. Using (C1), we have for δ ∈ (0, 1), s ∈ ℝ and r ∈ [1 − δ, 1] that


where we used that ∣δlogδ∣ ≤ C for δ ∈ (0, 1). Consequently, we have that G(r, s) ≥ −Cs∣ for δ ∈ (0, 1), s ∈ ℝ and r ∈ [1 − δ, 1]. A similar assertion holds for r ∈ [−1,−1 + δ]. Lastly, for ∣r∣ ≤ 1 − δ, we use (C2) to deduce that


and consequently G(r, s) ≥ −Cs∣ for all ∣r∣ ≤ 1 − δ and all s ∈ ℝ. Therefore, for all δ > 0, s ∈ ℝ and r ∈ ℝ, it holds that G(r, s) ≥ − Cs∣. Next, we define


By continuity of H(r) and (3.10), we can find a constant δ0 ∈ (0, r0 − 1) such that

H(r)>0 for r(1δ0,1+δ0),H(r)<0 for r(1δ0,1+δ0).

Then, it is clear that for ∣r∣ ≥ 2 r0, H(r) = 0 thanks to (3.10) satisfied by the extensions of fv and fφ. Meanwhile, for any δ ∈ (0, δ0), we see that if ∣r∣ ≤ 1 − δ0 < 1 − δ, then


which implies that

H(r)βlog,δ(r)C(1+|r|) for |r|1δ0.

On the other hand, as r0 > 1, for r ∈ [−2r0, −1 + δ0] ∪ [1 − δ0, 2r0], we use that βlog,δ(r) and H(r) have the same sign, so that their product H(r) βlog,δ(r) is non-negative. Hence, combining with the above analysis for the function G, we obtain the assertion (3.13).□

A. 4 Proof of Proposition 3.3


From the definition of β̂log,δ in (3.8), we infer that for r ≥ 1 − δ, δ ∈ (0, 1),


for some positive constant c independent of δ ∈ (0, 1). In a similar fashion, using βlog(δ − 1) < 0 and βlog(δ − 1) > 0, we have for r ≤ −1 + δ,


and when combined this yields


For the remaining case ∣r∣ ≤ 1 − δ, we employ the fact that βlog,δ(r) = βlog(r) and


to infer the existence of a constant c > 0 independent of δ ∈ (0, 1) such that

βlog(r)(r1)c for 0r<1,βlog(r)(r+1)c for 1<r0.

Hence, for δ ∈ (0, 1),


This completes the proof.□

B Well-posedness of the Brinkman system (2.9)

B. 1 Weak Solvability

Theorem 4

Let Ω ⊂ ℝd, d = 2, 3, be a bounded domain with C2,1-boundary Σ. Let cW1,r with r > d be given, fix exponent q such that q ∈ (1, 2] for d = 2 and q (65,2] for d = 3, and suppose η(⋅) and λ(⋅) satisfy (A3). Then, for any (f, g, h) ∈ Lq × L2 × (H1/2(Σ))*, there exists a unique weak solution (v, p) ∈ H1 × L2 to

div(2η(c)Dv+λ(c)div(v)IpI)+νv=finΩ, (B.1a)
div(v)=ginΩ, (B.1b)
Tc(v,p)n:=(2η(c)Dv+λ(c)div(v)IpI)n=honΣ, (B.1c)

in the sense


for all ΦH1. Furthermore, it holds that

vH1+pL2C(fLq+gL2+h(H1/2(Σ))) (B.2)

for a constant C > 0 depending only on Ω, q, η0, η1, λ0 and ν.


For (x1, …, xd)Ω, we define


Then, a short calculation shows that div(v0) = gΩ and Tc(v0, p0) = 0. Next, we will show there exist unique weak solutions (w, π) and (y, θ) to the systems

(P1)div(Tc(w,π))+νw=0 in Ω,div(w)=ggΩ in Ω,w=0 on Σ,(P2)div(Tc(y,θ))+νy=fνv0 in Ω,div(y)=0 in Ω,Tc(y,θ)n=hTc(w,π)n on Σ.

Then, one can check that the pair v := w + y + v0 and p := θ + π + p0 is the unique weak solution to (B.1).

B.1.1 Solvability of (P1)

By Lemma 2.1, there exists u := 𝓓(ggΩ) ∈ H01 satisfying


Let := −ν u + div(2 η(c) Du + λ(c) div(u)I) ∈ (H01), and consider the function space W0 := {f H01 : div(f) = 0 a.e. in Ω}. By the Lax–Milgram theorem, there exists a unique solution ŵW0 to


and satisfies

w^H1Cf~(H01)CuH1CgL2. (B.3)

Then, the function w := ŵ + u H01 satisfies div(w) = ggΩ a.e. in Ω and


i.e., w is the first component of the solution to (P1). The recovery of the unique pressure variable π L02 follows from the application of [51, p. 75, Lem. 2.2.2], which also yields

πL2CwH1CgL2. (B.4)

It is also clear that (w, π) constructed above is the unique solution to (P1), and for any ΦH1 it holds that


and so

Tc(w,π)n(H1/2(Σ))CgL2. (B.5)

B.1.2 Solvability of (P2)

We define the function space

Wdiv1,r:={fW1,r:div(f)=0 a.e. in Ω}.

By the Lax–Milgram theorem, there exists a unique solution y Wdiv1,2 to

F(Φ):=Ω2η(c)Dy:DΦ+νyΦfΦdx+hTc(w,π)n,ΦH1/2(Σ)=0 (B.6)

holding for all Φ Wdiv1,2 and satisfies

yH1C(fLq+h(H1/2(Σ))+gL2). (B.7)

By [51, p. 75, Lem. 2.2.2], there exists a unique pressure θ̂ L02 such that − ∇ θ̂ = F in the sense of distribution, with

θ^L2CF(H1)C(fLq+h(H1/2(Σ))+gL2). (B.8)

It remains to adjust this pressure by a uniquely defined constant c0 so that y and θ := θ̂ + c0 satisfy the boundary condition Tc(y, θ) n = hTc(w, π)n. Let q′ = qq1 denote the conjugate of q, where by the hypothesis it holds that q′ ≥ 2. From the distributional equation − ∇ θ̂ = F, we find that the weak divergence of 2 η(c) Dyθ̂ I satisfies

div(2η(c)Dyθ^I)=fνyLq, (B.9)

and so 2 η(c) Dyθ̂ I Ldivq. By [28, Thm. III.2.2], it holds that (2 η(c) Dyθ̂ I) n ∈ (W1/q,q(Σ))*, where by the generalised Gauss identity and (B.9) we have