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

We study a phase field model proposed recently in the context of tumour growth. The model couples a Cahn-Hilliard-Brinkman (CHB) system with a elliptic reaction-diffusion equation for a nutrient. The fluid velocity, governed by the Brinkman law, is not solenodial, 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 with singular potentials, specifically the double obstacle potential and the logarithmic potential, which ensures the phase field variable stay 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 may be of independent interests. As a consequence, included in our analysis is a weak existence result for a Darcy variant, and our work serve to generalise recent results on weak and stationary solutions to the Cahn-Hilliard inpainting model with singular potentials.


Introduction
Phase eld models [17,45] have recently emerged as a new mathematical tool for tumour growth, which o er 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 ngering 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].
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 uid velocity eld. The variable σ denotes the concentration of the nutrient and is assumed to evolve quasi-statically, while φ denotes the di erence in the volume fractions of the cells, with the region {φ = } representing the tumour cells and {φ = − } representing the host cells. The uid 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 e ect, 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 de nition of the stress tensor, the functions η and λ represent the shear and bulk viscosities, respectively. An interesting contrast with previous phase eld models in two-phase ows 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 nd the relation Ω Γv(φ, σ) dx = Ω div(v) dx = Σ v · n dH and the typical no-penetration boundary condition v·n = on Σ×( , T) will lead to a mean-zero compatibility condition for Γv(φ, σ), which can not be ful lled 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 ux ∇µ · n = φv · n and zero pressure p = on Σ T := Σ × ( , T) to circumvent the compatibility condition on Γv. For the Brinkman model, the compatibility issue can be avoided by prescribing the frictionless boundary condition Tn = on Σ T , as considered in [19][20][21][22]. Hence, we furnish (1.1) with the following initial and boundary conditions (1.2d) where ∂nf = ∇f · n is the normal derivative on Σ, φ 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 rst 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 rst derivative appears in (1.1d), is continuously di erentiable. The common example is the quartic potential ψ(s) = (s − ) . The purpose of the present paper is to provide a rst 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) = ( − r ) + I [− , ] (r) = ( − r ) for r ∈ [− , ], +∞ otherwise, and for the latter the logarithmic potential ψ log (r) = θ ( + r) log( + r) + ( − r) log( − r) + θc ( − r ) for r ∈ (− , ) with positive constants < θ < θc is the standard example. Let us brie y motivate the need to consider singular potentials such as ψ do and ψ log . The variable φ has a physical interpretation as the di erence between the volume fractions of the tumour cells and the host cells. As volume fractions are non-negative and bounded above by , φ should lie in the physical range [− , ]. This natural boundedness of φ becomes particularly important for physical models in uid ow and biological sciences, where the mass density of the mixture expressed as an a ne linear function of φ may become negative if φ strays out of [− , ]. A counterexample in [14] shows that for polynomial ψ, the phase eld variable φ can take values outside [− , ], 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 ow 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 L ( , T; H (Ω)). This is equivalent to controlling the mean value µ Ω := |Ω| Ω µ dx in L ( , 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 rst establishing the assertion φ Ω (t) ∈ (− , ) for all t > , which holds automatically for suitably chosen initial data thanks to the property of mass conservation φ Ω (t) = φ Ω ( ) for all t > .
When the Cahn-Hilliard component is a xed 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) ∈ (− , ), 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 = , σ = and fφ(φ) = λ(I − φ) for given non-negative λ ∈ L ∞ (Ω) and |I| ≤ 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 I [− , ] of ψ do , as well as showing that the source terms lead to the conclusion that φ Ω (t) ∈ (− , ) for all t > . 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 su ces to prescribe suitable conditions on bv, bφ, fv and fφ at φ = ± . Namely, we ask that bv and bφ are non-negative and vanish at ± , while fφ − fv is negative (resp. positive) at (resp. − ). 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 modi ed 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 di culties 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 v = − ν ∇p − (µ + χσ)∇φ 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 insu cient 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 rst 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 eld. This is interesting as the uid 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 eld 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 simpli ed model similar to the subsystem (1.1c)-(1.1e) with v = and χ = , 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φ) + ∆ φ − ∆ψ (φ) + χ∆σφ = Γφ (φ, σφ) in Ω, 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 eld, as regularity of the corresponding velocity eld is simply insu cient. 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 signi cant role in our work.

Main results . 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 = (a jk ) ≤j,k≤d , B = (b jk ) ≤j,k≤d ∈ R d×d , the scalar product A : B is de ned as the sum d j,k= a jk b jk . For ≤ p ≤ ∞, r ∈ ( , ∞), β ∈ ( , ) and k ∈ N, the standard Lebesgue and Sobolev spaces on Ω and on Σ are denoted by L p := L p (Ω), L p (Σ), W k,p := W k,p (Ω), and W β,r (Σ) with norms · L p , · L p (Σ) , · W k,p and · W β,r (Σ) , respectively. In the case p = we use the notation H k := W k, with the norm · H k . The space H is the completion of C ∞ (Ω) with respect to the H -norm, while the space H n is de ned as {w ∈ H : ∂nw = on Σ}. For the Bochner spaces, we use the notation L p (X) := L p ( , T; X) for a Banach space X with p ∈ [ , ∞], 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 and also the subspace of L -functions with zero mean value: The corresponding function spaces for vector-valued or tensor-valued functions are denoted in boldface, i.e., L p , W k,p , H k , L p (Σ) and W β,r (Σ). Furthermore, we introduce the space where div is the weak divergence. We now state three auxiliary lemmas that are used for the study of our model. The rst lemma concerns the solvability of the divergence problem, which will be useful for the mathematical treatment of equation (1.1a). Since we invoke Lemma 2.1 frequently with a = |Σ| Ω f dx n, we introduce where by the smoothness of the normal n, see (A ) 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.    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 su cient for the assertions of Lemma 2.2, as opposed to the Cregularity 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 B, and we mention that similar results can be found in [2] in the context of resolvent operators and H∞-calculus.
, be a bounded domain with C -boundary Σ and outer unit normal n. Let c ∈ W ,r with r > d be given and x exponent q ≤ r such that q > for d = and q ≥ for d = , and suppose the functions η(·) and λ(·) satisfy (A ). Then, for any f ∈ L q , g ∈ W ,q , h ∈ W − /q,q (Σ), there exists a unique with a constant C depending only on η , η , λ , ν, q, c W ,r and Ω.

. 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 .
Then, there exists a weak solution (φ, µ, σ, v, p) to (1.1)-(1.2) with the logarithmic potential ψ log in the sense of De nition 2.1 and satis es (2.14). Furthermore, for a.e. t ∈ ( , T), it holds that (2.15) where u = D(Γv(φ, σ)). . Let ρ ,ρ and u 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 u + u = . Then, setting φ := u − u , we consider the function where P, A > are the constant proliferation and apoptosis rates, so that proliferation occurs only at the interface region {− < φ < }. From the modelling, we have Γv = αΓ, Γφ = ρ S Γ, α := ρ − ρ , ρ S := ρ + ρ , and so we can set where bv(± ) = bφ(± ) = and It is also clear that bφ(s) and bv(s) satisfy (C ) and (C ). Furthermore, we point out that it is possible to have α < in the caseρ >ρ .

. The stationary problem
The

De nition 2.2. A quintuple
is a weak solution to the stationary CHB system with the double obstacle potential ψ do if the following properties hold: (e) equation (1.1a) holds a.e. in Ω, while (2.7a), (2.7d) and are satis ed for all Φ ∈ H and ζ ∈ H .
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 satis ed and In particular, the quintuple (φ, µ, σ, v, p) is a strong stationary solution.
We say that (φ, µ, σ, v, p) is a weak solution to the CHD system with the logarithmic potential ψ log if property (c ) in De nition 2.1 holds along with properties (g) and (h).

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.

. Approximation potentials and their properties . . Double obstacle potential
We point out that in order to use Lemma 2.2 the approximate potential should at least belong to C (R). Fix δ > , which serves as the regularisation parameter, and we de nê Formally, it is easy to see thatβ do,δ (r) → I [− , ] (r) pointwise as δ → , and so will serve as our approximation for the double obstacle potential. Let β do,δ (r) =β do,δ (r) = (r + ψ do,δ (r)) ∈ C (R) denote the derivative of the convex partβ do,δ , which has the form Then, it is clear that β do,δ is Lipschitz continuous with ≤ β do,δ (r) ≤ δ for all r ∈ 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 de ned as above. Then, there exist positive constants C and C such that for all r ∈ R and for all δ ∈ ( , / ), Aside from approximating the singular potential, it would be necessary to extend the source functions bv, bφ , fv and fφ from [− , ] to the whole real line. Since the solution variable φ is supported in [− , ] (see (c ) of De nition (2.1)), the particular form of extensions outside [− , ] does not play a signi cant role and we have the exibility 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 R such that fφ The latter implies fφ(r) − fv(r)r is strictly negative (resp. positive) for r > (resp. r < − ). For the functions stated in Example 2.1, we can consider following extensions: For r ∈ R, we set bv(r) = α max( , P( − r )) bφ(r) = ρ S max( , P( − r )), It is clear that (3.5) is ful lled. For α = , we see that fv(r) = and so fφ(r) − fv(r)r = fφ(r) satis es (3.6).
Hence, the extensions ful l (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φ modi ed as above. Uniform estimates will be derived in the next section and then we can pass to the limit δ → to infer the existence of a weak solution to (1.1) with the double obstacle potential in the sense of De nition 2.1.

. . Logarithmic potential
We de ne the convex part of ψ log and its corresponding derivative aŝ For xed δ > , the approximation of ψ log is the standard one: with convex part and corresponding derivativê 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 de ned as above. Then, there exist positive constants C , . . . , C , such that for all r ∈ R and for all < δ ≤ min( , θ/( θc)), Once again, we extend the source functions bv and bφ from [− , ] to the whole real line so that (3.5) is satis ed, and instead of (3.6) we consider extensions of fv and fφ from [− , ] to the whole real line that satisfy with r = . Then, using ρ S − α > it is clear that rfv(r) − fφ(r) ful ls (3.10).

. Existence of approximate solutions
To unify our analysis, we use the notation θ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 de nition of ψ δ (·), we can employ Lemma 2.2, and infer, for every δ ∈ ( , ), the existence of a weak solution quintuple for a.e. t ∈ ( , T) and for all Φ ∈ H and ζ ∈ H .

. 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.

. . Nutrient estimates.
Choosing ζ = σ δ in (3.12c) and using the non-negativity of h(·) leads to from which we deduce that σ δ is uniformly bounded in L ∞ ( , T; H ). Elliptic regularity additionally yields By a standard argument with the comparison principle, testing with ζ = −(σ δ )− and ζ = (σ δ − )+ will yield the pointwise a.e. boundedness of σ δ : Then, the continuous embedding H → L ∞ and the assumptions on bv, bφ, fv and fφ lead to . . Estimates from energy identity.
Thanks to (3.16), the function u := D(Γv) ∈ H satis es the estimate 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 ∞ ( , T; W ,p ) for any p ∈ ( , ∞). Henceforth, we drop the index δ and reuse the variable u.

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

Proof of Theorem 2 -Stationary solutions
As with the time-dependent case, we extend bv, bφ , fv and fφ from [− , ] to R such that fφ

. Approximation scheme
We consider a smooth functionĝ : R → [ , ] such thatĝ(r) = for r ≥ andĝ(r) = for r ≤ , and de ne F : L (Ω) → R as where C F is a positive constant to be speci ed later. Furthermore, we denote by γ(r, s) the function γ(r, s) := rΓv(r, s) − Γφ(r, s), and introduce, for δ ∈ ( , ), the cut-o operator T δ ∈ C , (R) de ned as It is clear that T δ and its derivative T δ are bounded independently of δ. Then, we seek a solution φ δ ∈ W := H n to the approximate system in Ω, where σ δ ∈ H is the unique non-negative and bounded solution to the nutrient subsystem and v δ is the rst component of the unique solution (v δ , p δ ) to the Brinkman subsystem in Ω, (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 δ ∈ ( , ). Then, we derive enough uniform estimates to pass to the limit δ → in order to prove Theorem 2. The new element in the analysis is how we treat the convection term.

. Preparatory results
For xed u ∈ W := H n it is clear that there is a unique solution σu ∈ H to the nutrient subsystem (4.3), and for any pair u , u ∈ W with corresponding solutions σu , σu . For the Brinkman subsystem (4.4), invoking Lemma 2.3 with c = u ∈ W, g = Γv(u, σu) ∈ H , f = (ψ δ (u) − ∆u)∇T δ (u) ∈ L , and h = , we can immediately deduce for any u ∈ W and δ ∈ ( , ) the existence of a unique strong solution (vu , pu) ∈ W , × W , . Although this is su cient 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 modi cation of the cut-o operator T δ . Proposition 4.1. Suppose T δ ∈ C , (R), and let (v , p ), (v , p ) ∈ W , × W , be the solutions to (4.4) corresponding to data u , u ∈ W, respectively. Then, there exists a positive constant C independent of the One example of the modi ed cut-o operator T δ ∈ C , (R) is the function Proof. For convenience we writeη = η(u ) − η(u ),ψ δ = ψ δ (u ) − ψ δ (u ),T δ = T δ (u ) − T δ (u ) andΓv := Γv(u , σu ) − Γv(u , σu ). Letŵ := D(Γv) ∈ H where by (4.5) it holds that Then, testing the di erence of (4.4) withv −ŵ yields The left-hand side can be bounded below by while the right-hand side can be bounded above by In light of the regularities u , u ∈ W and v , v ∈ W , , the polynomial growth of ψ δ leading to the following di erence estimate the Lipschitz continuity and boundedness of T δ , as well as (4.6), we nd that for some small constant ε > . Invoking Korn's inequality and (4.6) in (4.7) and combining with the above estimate gives For the pressure, we test the di erence of (4.4) withq := D(p) ∈ H which satis es and in turn we obtain This completes the proof.
The unique solvability of the nutrient and Brinkman subsystems in turn yields the following result. In particular, this shows that vn H + pn L ≤ C + ∇un L ψ δ (un) L + ∆un L ∇un L ≤ C δ , (4.8) for a positive constant C δ depending on δ but independent of n, thanks to the boundedness of {un} n∈N in W. Hence, for xed δ ∈ ( , ), there exist functions σu ∈ H , vu ∈ H and pu ∈ L , such that along a non-relabelled subsequence, 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 L . Together with ∆un ∆u in L and ∇un → ∇u in L , we obtain 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) =: Γv n yields

. Existence of approximate solutions
In this section we x δ ∈ ( , ), and de ne operators A , A : W → W * by Since β δ is bounded and β δ has sublinear growth, we deduce that the operator A 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 A is strongly continuous. Then, by [52,Thm. 27.6] the sum A = A + A is a pseudomonotone operator.
We now claim that A is additionally coercive over W, i.e., Let us rst treat the velocity term in Au, u W . By the de nition of T δ in (4.1), we see that uT δ (u)∇u = T δ (u)T δ (u)∇u = T δ (u)∇T δ (u) a.e. in Ω. Hence, by the trace theorem we have (4.9) Let uu := D (Γv u (u, σu)) ∈ H which satis es where the boundedness of Γv(u, σu) can be deduced using similar arguments in Section 3.3. Testing (4.4) for (vu , pu) with vu − uu 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 u ∈ W we obtain the following useful inequalities: 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, for the obstacle potential, since β do,δ (s) = , β log,δ (s) = β log (s) for s ∈ [− + δ, − δ]. Hence, as ≤ T δ (s) ≤ , we see that where for the second last inequality we used L'Hopital's rule to con rm that β log (s) (β log (s)) / is continuous in [− , ], hence bounded. Returning to (4.9), we deduce that for a positive constant C independent of u and δ. Now, computing Au, u W gives for a positive constant C independent of u and δ. Recalling that F(u) = C F for u L > |Ω|, and so, in choosing C F = C, we arrive at Invoking [52,Thm. 27.A], we deduce for every δ ∈ ( , ) the existence of a weak solution φ δ ∈ W to (4.2). Setting we see that the equation with right-hand side Thanks to the regularity φ δ ∈ W, v δ ∈ H , σ δ ∈ H , and the sublinear growth of β δ , we easily infer that f δ ∈ L . On the other hand, choosing ζ = in (4.15) implies that f ∈ L . Then, by arguing as in [33,Sec. 3.1], we obtain that µ δ ∈ W for all δ ∈ ( , ).

. Passing to the limit
In (4.16a) we take ζ ∈ H and apply integration by parts to get Passing to the limit δ → then yields (2.17). Meanwhile, (2.11) or (2.12) can be recovered in the limit δ → 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 δ ) δ∈( ,δ ) together with v constitutes a solution to the corresponding Brinkman system in the sense of (2.7a). From the de nition (4.14) of µ δ we observe from (4.4) that (v δ , p δ ) satis es for Φ ∈ H . For the last term, after integrating by parts and using T δ (φ δ )Φ → φΦ in L (Ω) and in L (Σ), µ δ → µ in L (Ω) and in L (Σ), we see that for all Φ ∈ H . Hence, passing to the limit δ → in (4.25) allows us to recover (2.7a) and thus the quintuple (φ, µ, σ, v, p) is a stationary solution in the sense of De nition 2.2.
Moreover, from the above estimates and weak lower semicontinuity of norms, we know that Then, from (2.17) and elliptic regularity, we deduce that In light of this improved regularity and the Sobolev embedding H ⊂ L ∞ , it is easy to see that (µ+χσ)∇φ ∈ L q where q < ∞ for d = and q = for d = . Invoking Lemma 2.3, we then infer v W ,q + p W ,q ≤ C, (4.28) which completes the proof.

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 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 re exive weak compactness arguments and [50,Sec. 8,Cor. 4], for δ → along a non-relabelled subsequence, it holds that for any r ∈ [ , ), The identi cation θ = div(φv) follows analogously as in Section 3.4, where the assertion (3.35) now holds for arbitrary λ ∈ L ( , T; L ) by the strong convergence ∇φ δ → ∇φ in L ( , T; L ) and the weak convergence v δ v in L ( , T; L ). In order to obtain uniform estimates for the chemical potential µ δ in L ( , T; L ), we again follow the argument in Section 3.4. Namely, we pass to the limit δ → in (2.7b) to obtain (3.36), and use the uniform boundedness of ψ δ (φ δ ) in L ( , T; L ) from (5.1) to obtain that the limit φ satis es the pointwise bound (3.38). Choosing ζ = in (3.36) leads to (3.39) and obtain by contradiction argument that (φ(t)) Ω ∈ (− , ) for all t ∈ [ , T].

A Properties of approximation potentials A. Proof of Proposition 3.1
Proof. As ψ do,δ is bounded for |r| ≤ + δ, it su ces to show (3.4a) holds for |r| > + δ. By Young's inequality it is clear that for r > + δ with δ ∈ ( , / ), and a similar assertion holds also for r < − − δ. This establishes (3.4a).
Then, one can check that the pair v := w + y + v and p := θ + π + p is the unique weak solution to (B.1).
Then, the function w :=ŵ + u ∈ H satis es div(w) = g − g Ω a.e. in Ω and i.e., w is the rst component of the solution to (P ). The recovery of the unique pressure variable π ∈ L follows from the application of [51, p. 75, Lem. 2.2.2], which also yields It is also clear that (w, π) constructed above is the unique solution to (P ), and for any Φ ∈ H it holds that and so Tc(w, π)n (H / (Σ)) * ≤ C g L . (B.5) B. . Solvability of (P ).

B. Strong solvability
Let us rst state the following auxiliary result for the Brinkman system with constant viscosities.

Lemma B.1.
Let Ω ⊂ R d , d = , , be a bounded domain with C -boundary Σ and outer unit normal n. Let η, λ and ν be positive constants, and x an exponent q such that q > for d = and q ≥ for d = .
Proof. We consider the functions analogous to v , p , w, π, y and θ de ned in the proof of Theorem B.1 for the case η(c) = η and λ(c) = λ, whilst reusing the notation. It is clear that v W ,q + p W ,q ≤ C g W ,q .
By the assumption on the exponent q, it holds that W ,q ⊂ H and W ,q ⊂ L . Hence, the di erence w := w−z and π := π − ϕ constitutes a weak solution in H × L to the system Unique solvability implies ( w, π) = ( , ) and hence the unique weak solution (w, π) to (P ) with constant viscosities is in fact a strong solution satisfying w W ,q + π W ,q ≤ C f L q + g W ,q .

B. Proof of Lemma 2.3
Fix c ∈ W ,r with exponent r > d and data f ∈ L q , g ∈ W ,q and h ∈ W − /q,q (Σ), where the exponent q ≤ r satis es q > for d = and q ≥ for d = .
Note that for d = , we have the Sobolev embedding for all ψ ∈ H . Since q > and − q ∉ Z, surjectivity of the trace operator yields the existence of an extension h ∈ W ,q (Ω) of h ∈ W − /q,q (Σ) such that h W ,q ≤ C h W − /q,q (Σ) . Furthermore, as r > d, we have the Sobolev embedding W ,q ⊂ L qr r−q and c ∈ W ,r ⊂ C , − d r (Ω). Hence, it is easy to see that η(c) h ∈ W ,q and by the trace theorem η(c) h ∈ W − /q,q (Σ). Next, we de ne the exponent s = r +r < so that At this point the analysis is divided into two cases: B. . Case 1 (q ≤ s = r +r ).
We have q < and q −q ≤ r, and so , which implies k ∈ L q . From (B.16), we see that (v, η(c)p ) ∈ H × L is a weak solution to the following system with constant viscosity where k ∈ L q , g ∈ W ,q and η(c) h ∈ W − /q,q (Σ). By Lemma B.1, it holds that (v, η(c)p ) ∈ W ,q × W ,q is in fact a strong solution satisfying v W ,q + η(c)p W ,q ≤ C k L q + g W ,q + η(c) h W − q ,q (Σ) ≤ C + ∇c L r f L q + g W ,q + h W − q ,q (Σ) . By the assumption on η, we easily infer that p L q ≤ C η(c)p L q and p W ,q ≤ C η(c)p L q + η(c)∇( η(c)p ) L q + p η (c) η(c) ∇c L q ≤ C η(c)p W ,q + C ∇c L r p L qr r−q ≤ C η(c)p W ,q + ∇c L r p L ≤ C + ∇c L r f L q + g W ,q + h W − q ,q (Σ) .
where we used q −q ≤ r to deduce that qr r−q ≤ . De ning p :=p + λ(c)g and making use of the fact div(v) = g we see from (B.15) that (v, p) satis es the Brinkman system (2.9). Furthermore, from the estimate p W ,q ≤ p W ,q + λ(c)g W ,q ≤ p W ,q + C ∇c L r g L qr r−q + C ∇g L q ≤ p W ,q + C + ∇c L r g W ,q where we used r > d to deduce that − d q ≥ − d(r−q) qr and leads to the Sobolev embedding W ,q ⊂ L qr r−q , we nd that (v, p) also satis es the estimate (2.10).
B. . . Case 2 (q > s = r +r ). In this case, using s < and r = s −s , we see that .
Let t be the exponent de ned by then by Sobolev embedding it holds that (v,p) ∈ W ,t × L t . Since r > d, it follows that t > . The idea is to use this improved regularity as a starting point of a bootstrapping argument to show k ∈ L q . Once we have k ∈ L q , following the argument in Case 1 then gives the desired assertion. Let s < t be the exponent de ned by Since r > d, we have s > s, and from (B.16), k L min(s ,q) ≤ C f L q + C v L t + C Dv L t + p L t ∇c L r .
If q ≤ s then k ∈ L q and the proof is complete, otherwise if s < q, we deduce instead that (v,p) ∈ W ,s × W ,s ⊂ W ,t × L t where