Cahn–Hilliard–Brinkman model for tumor growth with possibly singular potentials

We analyze a phase field model for tumor growth consisting of a Cahn–Hilliard–Brinkman system, ruling the evolution of the tumor mass, coupled with an advection-reaction-diffusion equation for a chemical species acting as a nutrient. The main novelty of the paper concerns the discussion of the existence of weak solutions to the system covering all the meaningful cases for the nonlinear potentials; in particular, the typical choices given by the regular, the logarithmic, and the double obstacle potentials are admitted in our treatise. Compared to previous results related to similar models, we suggest, instead of the classical no-flux condition, a Dirichlet boundary condition for the chemical potential appearing in the Cahn–Hilliard-type equation. Besides, abstract growth conditions for the source terms that may depend on the solution variables are postulated.


Introduction
Cancer is nowadays still one of the main diseases causing death worldwide.Beyond doubt, the understanding of the development of solid tumor growth is one of the major challenges scientists have to face in the current century.Moreover, it is now, more than ever, apparent that only interdisciplinary efforts may enable us to gain deeper insights into cancer development mechanisms.In this scenario, mathematics could play a crucial role, since multiscale mathematical modeling provides a quantitative tool that may help in diagnostic and prognostic applications (see, e.g., the seminal book [8]).Among others, mathematics has two decisive advantages: the first one is that of being able to select particular mechanisms that may be more relevant than others, while the second one is that of being able to foresee and make predictions that may be precious for medical practitioners, without inflicting any harm to the patients.Furthermore, the extremely fast development of computational methods for the solution of nonlinear PDEs opens the doors for a direct interaction between the experimental methods used by physicians and the more theoretical mathematical ones: indeed, advanced numerical solvers may be implemented as a supporting tool in clinical therapies.Recently, lots of phase field models modeling tumor growth have been proposed: a brief description of the state of the art will be provided later on.As biological materials like tumor agglomerates exhibit viscoelastic properties, we prescribe a velocity equation of Brinkman type.
Let Ω ⊂ R 3 be a spatial domain in which the tumor is located, T > 0 be a fixed final time, and set Q := Ω × (0, T ), and Σ := ∂Ω × (0, T ).Then the system under investigation in this paper is a Cahn-Hilliard-Brinkman model related to tumor growth and reads as follows: where the coefficients m and n are positive functions and the viscous Brinkman stress tensor T is defined as Here, the standard notation is used for the symmetrized gradient of the velocity field v which represents the volumeaveraged velocity field of the mixture with permeability ν.Moreover, p is the pressure, η and λ are nonnegative constants denoting the shear viscosity and the bulk viscosity, respectively, and I ∈ R 3×3 is the identity matrix.Although several choices are possible, we endow the above system with the following boundary and initial conditions: on Σ, (1.9) ϕ(0) = ϕ 0 , σ(0) = σ 0 in Ω, (1.10) where n denotes the outer unit normal vector to ∂Ω, and ∂ n the associated outward normal derivative.The other variables of the system are ϕ, µ, and σ.In relevant cases, the phase variable ϕ is an order parameter taking values in [−1, 1] that represents the difference between the volume fractions of tumor cells and healthy cells.It allows us to keep track of the evolution of the boundary of the tumor, since the level sets {ϕ = 1} := {x ∈ Ω : ϕ(x) = 1} and {ϕ = −1} describe the region of pure phases: the tumorous phase and the healthy phase, respectively.The second variable µ denotes the chemical potential related to ϕ as in the framework of the Cahn-Hilliard equation.We postulate the growth and proliferation of the tumor to be driven by the absorption and consumption of some nutrient σ (usually oxygen).Whenever 0 ≤ σ ≤ 1, σ ≃ 1 represents a rich nutrient concentration, whereas σ ≃ 0 a poor one.The functions m(ϕ) and n(ϕ) are nonnegative mobility functions related to the phase and to the nutrient variables, respectively.The positive physical constants ǫ and ν are connected to the interfacial thickness and surface tension, while the nonnegative constant χ stands for the chemotactic sensitivity.Finally, S ϕ and S σ denote nonlinearities representing some source terms that account for the mutual interplay between tumor, healthy cells, and nutrients.For further details concerning the modeling, we refer to [16] and the references therein.
Concerning the boundary conditions, we point out that (1.8) can be understood as a no-friction boundary condition for the velocity field v.This is a common request for similar systems, see, e.g., [9,10,23], as it does not enforce any compatibility condition on the velocity source term g in (1.2), which is otherwise needed if one assumes a no-slip boundary condition like v = 0 on Σ or the no-penetration boundary condition v • n = 0 on Σ.Both these conditions entail that Ω g = 0, which is not ideal from the modeling perspective.Let us also notice that the vector field v is not solenoidal (as typically in fluid-type problems), which entails challenges from the mathematical viewpoint.
As a common denominator of a more general Cahn-Hilliard equation, in (1.4) F ′ represents the (generalized) derivative of a double-well shaped nonlinearity F .Prototypical examples are the regular, the logarithmic, and the double obstacle potentials.These read, in the order, as ) for some positive constants θ < θ 0 and c.Besides, in the case of nonregular potentials like the double obstacle (1.13), the second equation (1.4) has to be read as a differential inclusion.
The main novelty of this paper is the prescription of the Dirichlet boundary condition µ = µ Σ for the chemical potential on Σ, in contrast to the standard homogeneous Neumann (no-flux) boundary condition ∂ n µ = 0 on Σ (see, e.g., [9,10,23]).The limitation behind the latter choice regards the nonlinear potentials that can be considered: in all of the aforementioned papers, the authors were forced to restrict the analysis to regular potentials, possibly of just quadratic growth at infinity; singular potentials like (1.12) or (1.13) were excluded from the analysis.These, however, are actually more relevant physically, since, if solutions to the system exist, then the condition ϕ ∈ [−1, 1] is automatically fulfilled.The restriction of the admitted potentials originates from the presence of the source term S ϕ in the Cahn-Hilliard equation (1.1).Roughly speaking, for proving that µ ∈ L 2 (0, T ; H 1 (Ω)), the energetic approach provides just a control on ∇µ in L 2 (Q).The classical approach then requires the employment of the Poincaré-Wirtinger inequality along with a control of the spatial mean of µ in L 2 (0, T ).This can be achieved by comparison in equation (1.4), provided that F is regular and its derivative F ′ possesses a prescribed growth, which leads to the choice of potentials of polynomial type.Therefore, the novelty of our work revolves around the different boundary condition, which allows us to apply another version of Poincaré's inequality to establish immediate control of µ ∈ L 2 (0, T ; H 1 (Ω)) from the bound ∇µ ∈ L 2 (Q) without the need of an additional control of the spatial average of µ; in this way, unpleasant growth restrictions for the potential can be avoided (see also [14]).
Without claiming to be exhaustive, let us now review some literature connected with the system (1.1)-(1.10).It is well known that the Brinkman law interpolates between the Stokes and Darcy paradigms, and it has become rather popular in recent times, see [9,10,23].For tumor growth models both descriptions seem reasonable, because the associated Reynolds number is very small.Formally, we recover the Darcy limit when η ≡ λ ≡ 0, and ν > 0, where the boundary condition (1.8) yields that p = 0 on Σ; similarly, the Stokes limit is obtained when η, λ > 0, and ν = 0.The Stokes equation was suggested, e.g., in [5,11], by approximating the tumor as a viscous fluid, while Darcy's law describes a viscous fluid permeating a porous medium represented by the extracellular matrix and accounts for the inclination of cells to move away from regions of high compression, see, e.g., [15,16].Stationary approximations for system (1.1)-(1.10)are popular as well, and we mention [13,14,18] and the references therein, where just polynomial-type potentials were considered.To cope with the case of singular potential, some authors (see [6]) suggested to include suitable relaxations.Besides, we refer to [19,20,28] and the references therein, for related nonlocal versions, to [1,17,24] for the additional coupling with elasticity, and to [27] for the coupling with the Keller-Segel equation.

Biological Examples and Modeling Considerations
Before diving into the mathematical details, let us outline some physically meaningful choices for the source terms S ϕ and S σ introduced above.Linear kinetic.A first typical form for S ϕ and S σ was proposed by Cristini et al. in [8] motivated by linear kinetic: (1.14) In the above expressions, P, A, B, and C denote nonnegative constants related to biological quantities that are, in the order, the proliferation rate of tumoral cells by consumption of nutrient, the apoptosis rate, the consumption rate of the nutrient with respect to a preexisting concentration σ c , and the nutrient consumption rate.As for the function , it denotes an interpolation function between −1 and 1 with the property that (−1) = 0 and (1) = 1.Roughly speaking, weights the corresponding mechanism compared to the amount of cancer located in that region and "turns off" the associated mechanism when the tumor is not present.
Linear phenomenological laws for chemical reactions.Another approach was proposed by accounting for linear phenomenological laws for chemical reactions by A.
Hawkins-Daarud et al. in [22], where the following form was suggested: where P stands for a suitable nonnegative proliferation function.
In the forthcoming analysis, we will proceed in an abstract fashion without prescribing explicit structures for the source terms, but just postulating suitable growth conditions.Those are straightforwardly fulfilled by the above cases but this last one.Namely, we cannot allow a linear growth of the sources with respect to µ.Nevertheless, let us claim that the special form (1.15) can still be considered provided to adjust some estimates accordingly (cf.[6,7]).
The above system (1.1)-(1.10)(see, e.g., [16]) is naturally associated to a free energy E with following the structure: with N specified by (1.17) Here, the first two terms in E yield the well-known Ginzburg-Landau energy modeling phase segregation and adhesion effects, while N stands for the chemical free energy density, respectively.For convenience, let us immediately set a specific notation to denote the partial derivatives of the free energy N with respect to the variables ϕ and σ: as a direct computation shows.This notation will turn very convenient for the estimation procedure later on.
As they will play any role from the mathematical viewpoint, we will set m(•) ≡ n(•) ≡ 1, and ǫ = 1 in our discussion.However, let us claim that our well-posedness result might be proven also for suitable nonconstant, albeit non-degenerate, mobilities (see [9]).

Notation, Assumptions and Main Results
Throughout the paper, Ω is a bounded and connected open subset of R 3 (the twodimensional case can be treated in the same way) having a smooth boundary Γ := ∂Ω.In the following, |Ω| denotes the Lebesgue measure of Ω.Similarly, we write |Γ| for the two-dimensional Hausdorff measure of Γ.Given a final time T > 0, we set, for every t ∈ (0, T ], (2.1) Given a Banach space X, we denote by • X , X * , and •, • X , its norm, its dual space, and the associated duality pairing, respectively.As for the notation of norms, some exceptions will be utilized in the sequel.Moreover, the symbol used for the norm in some space X is adopted for the one in any power of X.Another standard notation that we employ concerns vectors, or vector-valued functions, which are denoted by bold symbols; for instance, 0 stands for the zero vector in R 3 , and v = (v 1 , v 2 , v 3 ) ∈ R 3 .For 1 ≤ q ≤ ∞ and k ≥ 0, we indicate the usual Lebesgue and Sobolev spaces on Ω by L q (Ω) and W k,q (Ω), with the standard abbreviation H k (Ω) := W k,2 (Ω).The norm in L q (Ω) is simply denoted by • q , and the same symbol is used for the norms in the analogous spaces constructed on Q, Γ and Σ, if no confusion can arise.Then, we introduce the shorthands , and W := {v ∈ H 2 (Ω) : ∂ n v = 0 a.e. on Γ}, and endow these spaces with their natural norms.For simplicity, we write • instead of • H . Besides, the space H will be identified with its dual, so that we have the following continuous, dense, and compact embeddings: ) is a Hilbert triplet that will be used as well.We observe at once the compatible embeddings and we recall that Finally, the same symbols written with boldface characters denote the corresponding spaces of vector-valued functions.So, we have for instance that L p (Ω) = (L p (Ω)) 3 , H = H 3 and V = V 3 .However, since H and V are powers of H and V , we simply write • and • V instead of • H and • V , respectively.Moreover, for given matrices A, B ∈ R 3×3 , we define the scalar product The following structural assumptions will be in order in our analysis.The double-well potential F introduced in the examples (1.11)-(1.13) is replaced by a more general one.Indeed, we can make the following assumptions: (2.5) For the source terms we assume: for some constant Θ > 0 and every r, s ∈ R . (2.8) Finally, the permeability, viscosity and sensitivity constants ν, η, λ and χ are requested to satisfy ν, η ∈ (0, +∞) and λ, χ ∈ [0, +∞). (2.10) It is well known that β is a maximal monotone graph in R ×R with corresponding domain D(β), and that 0 ∈ β(0).It is worth noticing that in the case β ∈ C 1 (R) it follows that it is single-valued, and we can write β = β ′ as well as F ′ = β + π.Here, we immediately observe that all of the standard potentials (1.11)-(1.13)fulfill (2.4)-(2.6),as well as that the biologically relevant examples given above for the source terms satisfy (2.7)-(2.9).
(2.13)However, it is convenient to transform the Dirichlet inhomogeneous boundary condition for µ into a homogeneous one by performing a change of variable for the chemical potential.
To this end, we introduce the harmonic extension of the boundary datum µ Σ , i.e., the function h : Q → R defined as follows: We notice at once that (2.12) ensures that h enjoys at least the regularity Thus, upon setting µ := µ − h, the boundary condition for µ in (1.9) now becomes the homogeneous Dirichlet condition µ = 0 on Σ.However, for the sake of simplicity, we proceed with abuse of notation and still denote by µ the above difference µ between the chemical potential and h.This change of notation obviously slightly modifies the equations: the right-hand side of (1.1) has to be adapted, and we have to rewrite (1.4) as where the differential inclusion arises from the possible multi-valued nature of the nonlinearity F , in accordance with (2.4)- (2.6).For this reason, and in order to clarify the meaning of the equation, we state the new problem to be dealt with in a precise form.
In particular, the equations (1.1), (1.3) and (1.5), and the corresponding boundary conditions, are replaced by variational equations (also owing to the Leibniz rule for the divergence and (1.2)), and the homogeneous Dirichlet condition for µ is enforced by the forthcoming (2.18).
By recalling (1.6)-(1.7)and (1.17)-(1.18),and setting m(•) ≡ n(•) ≡ 1 and ǫ = 1 as announced in the Introduction, we look for a sextuple (v, p, µ, ϕ, ξ, σ) enjoying the regularity properties ) that satisfies for every φ ∈ V 0 and a.e. in (0, T ) , (2.24) for every z ∈ V and a.e. in (0, T ) , (2.25) for every ζ ∈ V and a.e. in (0, T ) , (2.27) as well as the initial conditions ϕ(0) = ϕ 0 and σ(0) = σ 0 . (2.28) Here is our main result: As further regularity is concerned, the components ϕ and ξ of every solution satisfy (in the two-dimensional case the summability exponent 6 can be replaced by any q ≥ 1), as shown in the forthcoming Remark 4.1.Finally, we notice that, in view of the very low regularity at disposal for the nutrient variable σ and the velocity field v, the uniqueness of weak solutions is not to be expected.
We continue this section by listing some tools that will be useful later on.We first recall Young's inequality where q ′ denotes the conjugate exponent of q given by the identity (1/q) + (1/q ′ ) = 1.We repeatedly use it, mainly with q = q ′ = 2.We also account for Hölder's inequality, as well as for the following Sobolev, compactness, Poincaré, and Korn inequalities: Here, the constants C S , C P , and C K , depend only on Ω, while C δ depends on δ, in addition.In (2.35), the notation (1.7) is used.
Next, we present three auxiliary results that will be used in the sequel.The first one is stated in a more general setting as an exercise in [12, Ex.III.3.5], but it readily follows as a corollary from [12,Thm. III.3.1].The second one is related to the Stokes resolvent operator and is a particular case of [2, Thm.3] (which is an extension of the results in [21]).Finally, the last one regards the trace operator.
Lemma 2.3.There exists a constant C that depends only on Ω such that for every f and a satisfying f ∈ H, a ∈ H 1/2 (Γ), and there exists some u ∈ V satisfying (2.37) Lemma 2.4.Assume that f ∈ H and f ∈ V .Then, there exists a unique pair (v, p) Moreover, the mapping Lemma 2.5.The trace operator maps L ∞ (0, T ; H)∩L 2 (0, T ; V ) into L 4 (0, T ; L 2 (Γ)), and it holds the estimate , where the constant C depends only on Ω, and where v also denotes the trace of v on Γ.
Proof of Lemma 2.5.As we did not find a precise reference, we provide a sketch of the proof.We denote by C 1 , C 2 ,. . .constants that depend only on Ω.We introduce the real interpolation space (by the way, the Besov space (2.39) The trace operator v → v| Γ maps B into L 2 (Γ) and is linear and continuous.In the half-space case Ω = R 3 + , this can be deduced, e.g., from formula (I.17) (with the notation Y (1, R 3 + ) for (2.39)) in the paper [30], where it is also shown that the operator maps B onto L 2 (Γ).As usual, the result is then extended to the general case by using local charts and a partition of unity.This leads to the estimate On the other hand, the interpolation inequality L ∞ (0,T ;H) , and (2.38) directly follows by taking the 4th powers and integrating over (0, T ).
Remark 2.6.We notice that the application of the trace estimate in [12, Thm.II.4.1] (whose proof is left to the reader), with the parameters therein being chosen as r = q = 2, m = 1, n = d = 3, λ = 0, produces a similar trace inequality, namely, Besides, let us state a general rule concerning the constants that appear in the estimates to be performed in the following.The small-case symbol c stands for a generic constant whose actual value may change from line to line, and even within the same line, and depends only on Ω, the shape of the nonlinearities, and the constants and the norms of the functions involved in the assumptions of the statements.In particular, the values of c do not depend on the parameters ε and k that will be introduced in the next section.A small-case symbol with a subscript like c δ (specifically, with δ = ε) indicates that the constant may depend on the parameter δ, in addition.On the contrary, we mark precise constants that we can refer to by using different symbols (see, e.g., (2.32) and (2.37)).
The next sections aim to rigorously prove Theorem 2.1.A standard approach to guarantee the existence of solutions to similar Cahn-Hilliard type systems is based on suitable approximation procedures that can be schematized as follows: first, one has to regularize the possibly singular nonlinearity F , with the classical choice being the Yosida regularization that depends on a parameter, say, ε > 0.Then, for every fixed ε > 0, one further discretizes the system in space via the Galerkin method and solves a family of finite-dimensional problems that depend on a parameter k ∈ N. Next, one has to provide some rigorous estimates, which are independent of both ε and k.From these estimates, using weak and weak star compactness arguments, one can find a suitable subsequence and eventually pass to the limit as k → ∞ and as ε → 0, thus showing that the obtained limits yield a solution to the original system.
The essence of the proof of Theorem 2.1 can be roughly schematized by the abovementioned steps, but we had to face a major obstacle arising from the different boundary conditions in (1.9).The latter prevent the solvability of the ODE system originating from the Galerkin scheme.In fact, due to (1.9), one is in the approximation naturally led to consider Schauder bases (cf.(3.18) and (3.19)) for the Laplace operator with homogeneous Neumann and Dirichlet boundary conditions for ϕ, σ and µ, respectively.Despite of being natural, this choice completely impedes the solvability of the discrete problem, because there occur inner products between the elements of the different two bases.To overcome this intrinsic difficulty, we introduce an intermediate approximation step, which consists in adding further regularizing terms at the level ε > 0 that somehow dominate the mixed terms and enable us to solve the Galerkin system (cf.(3.23)-(3.28)).Finally, we pass to the limit as ε ց 0 as anticipated above, thus proving the theorem.

Approximation
In this section, we introduce and solve a proper approximating problem depending on the parameter ε > 0. First of all, we replace the functional β and the graph β by their Moreau-Yosida regularizations β ε and β ε , respectively (see, e.g., [4, pp. 28 and 39]).Then, we set and recall that classical theory of convex analysis entails the following facts: for every M > 0 there exist C M > 0 and ε M > 0 such that The only nonobvious fact is the coercivity property in (3.4).In this direction, let us fix M > 0 and observe that our assumption (2.5) on β implies that It then follows that where s * is the minimum point, namely, s * = r/(1 + 4Mε).Hence, we have that (3.4) with F ε replaced by β ε (with an obvious choice of ε M ).Then, (3.4) itself follows, since π is Lipschitz continuous.
Besides this regularization, we replace g and h by smoother functions g ε and h ε satisfying g ε ∞ ≤ c , and g ε → g a.e. in Q as ε ց 0 , (3.5) For simplicity, we do not enter in the details concerning the construction of the regularizations above.However, let us mention that standard mollification arguments are enough to get the desired properties prescribed by (3.5) and (3.6), respectively.Moreover, as anticipated, we introduce artificial viscous terms in some of the equations.We prefer to present all of them in their variational form.The approximating problem thus consists in finding a quintuple (v ε , p ε , µ ε , ϕ ε , σ ε ) that satisfies the regularity properties and solves the variational equations for every ζ ∈ V and a.e. in (0, T ) , (3.12) for every φ ∈ V 0 and a.e. in (0, T ) , ε for every z ∈ V and a.e. in (0, T ) , for every ζ ∈ V and a.e. in (0, T ) , as well as the initial conditions Let us incidentally notice that the presence of a selection ξ (cf.(2.26)) is no longer needed in the approximating ε-problem due to the regularity at disposal for F ε defined by (3.1).Moreover, we can simply write Remark 3.1.We notice that, due to (3.10) and the initial condition for ϕ ε , the regularity for ∂ t µ ε given in (3.9) and the initial conditions for µ ε are equivalent to respectively.
Theorem 3.2.Under the assumptions of Theorem 2.1 and with the above notation, the approximating problem (3.12)-(3.17)has at least one solution The rest of this section is devoted to the proof of this theorem.The method we use starts from a discrete problem based on a Faedo-Galerkin scheme.To this end, we introduce the nondecreasing sequences {λ j } and {λ 0 j } of eigenvalues and the corresponding complete orthonormal sequences {e j } and {e 0 j } of eigenfunctions of the eigenvalue problems for the Laplace operator with homogeneous Neumann and Dirichlet boundary conditions, respectively.Namely, we have that − ∆e j = λ j e j in Ω and ∂ n e j = 0 on Γ , (3.18) − ∆e 0 j = λ 0 j e 0 j in Ω and e 0 j = 0 on Γ , for j = 1, 2, . . ., as well as the normalization conditions Ω e i e j = Ω e 0 i e 0 j = δ ij for every i and j , with the standard Kronecker symbol δ ij .Moreover, if we set, for k = 1, 2, . . ., then the unions of the these spaces are dense in V and V 0 , respectively, and both are dense in H as well.We notice at once that all of the above eigenfunctions are smooth since Ω is smooth, that λ 1 = 0, and that e 1 = |Ω| −1/2 .Then, the discrete problem related to k consists in finding a quintuple (v k , p k , µ k , ϕ k , σ k ) with the regularity specified by that solves the system ε ε ) for every ζ ∈ V , φ ∈ V 0 k , z, ζ ∈ V k , and t ∈ [0, T ), and that fulfills the initial conditions and Existence for the discrete problem.The first aim of ours is to show the existence of at least one solution (we do not care about uniqueness, since it is not needed).The method relies on a proper application of Lemma 2.4.Besides, let us point out that the idea of employing the Stokes resolvent to express the velocity field v k in terms of the other variables ϕ k , µ k , and σ k is largely inspired by [9] (see also [10,23]).For a while, the symbols ϕ k , µ k and σ k denote independent variables.To every triplet (ϕ k , µ k , σ k ) ∈ V k × V 0 k × V k , we associate the vector-valued function and we notice that f k only depends on time through the continuous function h ε , while the dependence on space occurs just through the eigenfunctions and their gradients, which are smooth.In particular, on the one hand, if we read f k as a function of the coefficients of ϕ k , µ k , σ k (with respect to the bases just chosen), and t, then we see that it is continuous.On the other hand, for every t ∈ [0, T ], we are allowed to apply Lemma 2.4 with f = f k (t) and f = g ε (t).Since the mapping Ψ is linear, continuous, and time independent, and since g ε and h ε are continuous, this yields a pair of functions that are continuous with respect to the coefficients of ϕ k , µ k , σ k , and t.This observation is made to ensure the continuity of the functions that rule the system of ODE's we are going to introduce.Now, we let ϕ k , µ k and σ k depend on time.To every triplet (ϕ we associate the function f k still given by (3.29) and, for every t ∈ [0, T ], we apply Lemma 2.4 as before.We obtain two functions, which we still term Ψ 1 (f k , g ε ) and Ψ 2 (f k , g ε ) with an abuse of notation, that, according to the lemma, belong to L ∞ (0, T ; H 2 (Ω)) and L ∞ (0, T ; V ), respectively.By construction, the pair of functions (v k , p k ) := (Ψ 1 (f k , g ε ), Ψ 2 (f k , g ε )) solves the equations (3.23)-(3.24)corresponding to the given triplet (ϕ k , µ k , σ k ).Therefore, the whole problem (3.23)-(3.28) is equivalent to the problem of finding a triplet (ϕ k , µ k , σ k ) with the regularity specified in (3.22) that solves for every φ ∈ V 0 k , z, ζ ∈ V k and t ∈ [0, T ), with f k given by (3.29), and satisfies the initial conditions in (3.28).We show that this problem has at least one solution.To this end, we represent the unknowns in terms of the bases of the spaces V k and V 0 k , i.e., which are the true unknowns.In terms of these coefficients, the discrete problem takes the form with some continuous functions , and the initial conditions for ϕ k , µ k and σ k are trivially derived from (3.28).By replacing ϕ ′ k in (3.33) using (3.34) (recall that now ε > 0 is fixed), we obtain a standard Cauchy problem for a 3k-dimensional nonlinear ODE system ruled by a continuous function.This allows us to apply the Cauchy-Peano theorem, which ensures the existence of at least one local solution.This local solution can be extended to a maximal solution, which provides a maximal solution (ϕ k , µ k , σ k ) to (3.23)-(3.28)defined in the interval [0, T k ) for some T k ∈ (0, T ].We claim that this solution is bounded (as required) and global, i.e., that T k = T .The proof relies on the estimate which we are going to prove in the next lines.Since (3.20) implies that and similarly for the other two components, (3.36) shows that the R 3k -valued function ( ϕ k , µ k , σ k ) is bounded.Then, maximality also implies that the solution is global.
Next, we perform a number of a priori estimates that will allow us to let k tend to infinity and to show that the approximating problem (3.12)-(3.17)has at least one solution.The first of these estimates proves the validity of (3.36), in particular.Remark 3.3.We note that the initial values ϕ k (0) and σ k (0) are the H-projections of ϕ 0 and σ 0 onto V k .To deal with them and for other purposes, it is worth noting a property of the H-projection v k of a generic element v ∈ V onto V k .We have that v k ≤ v , and we derive a similar inequality for the gradients.By definition, for j = 1, . . ., k we have that Ω v k e j = Ω ve j , whence, using (3.18), also

By linear combination, it follows that
and we conclude that the H-projection v k coincides with the V -projection of v.
By noting that the choice w = v k is admitted in the above identity, and collecting everything, we conclude that Therefore, in particular, we have the inequalities Now let q ∈ [1, +∞] and v ∈ L q (0, T ; V ), and define v k : Q → R as follows: for a.a.
Moreover, if q < +∞, then we also have that v k → v strongly in L q (0, T ; V ).Indeed, q for a.a.t ∈ (0, T ) , so that one can apply the Lebesgue dominated convergence theorem.Clearly, everything can be repeated for the space V 0 and the projection on V 0 k .
First a priori estimate.We recall that we do not yet know that T k = T .For every t ∈ (0, T k ), we apply Lemma 2.3 with the choices by observing that the assumptions (2.36) are satisfied, and we term u k (t) the function given by the lemma.Avoiding writing the time t for a while, and recalling (2.37) and (3.5), we have that Then, we test (3.23) by ζ = v k − u k , with the aim of removing the pressure from the first estimate.This idea has been introduced in [9], and it relies on the identities Dv k : Besides, by also integrating over (0, t) with respect to time, where t ∈ (0, T k ) is arbitrary, we obtain that 2η as well as (1.18), testing (3.27) by N σ (ϕ k , σ k ) and integrating with respect to time, we also obtain that Next, we observe that the Korn inequality (2.35) implies that 2η Finally, as for the term involving F ε , we apply (3.4) with M = 1, and obtain that provided that ε is small enough (as in the lemma): from now on, it is understood that ε satisfies this smallness condition.Regarding the right-hand side, we start by treating the nontrivial terms coming from the identity (3.40).The symbol δ denotes a positive parameter whose value is chosen later on.By recalling (3.39), we have that Next, with the help of the Hölder, Sobolev, Poincaré, and Young inequalities, we obtain that For the next term, we also apply the compactness inequality (2.33) to N σ (ϕ k , σ k ) and have that Similarly, using (3.6) and the Sobolev inequality (2.32) for h ε and v k , we have that We notice that the L 1 norm of the function s → h ε (s) 2 V is bounded by a constant independent of ε by (3.6).Now, among the volume integrals that come from (3.41)-(3.43)and should be estimated, just the last one on the right-hand side of (3.42) needs some treatment.Indeed, all the others can easily be dealt with by virtue of the Young inequality, possibly combined with other estimates, like (2.8) or (2.34), without any difficulty.We have that where we owe to (3.6) for the last inequality.Now, we move to the surface integrals.We have that Finally, all of the terms involving the initial values can easily be estimated by accounting for (3.38) and recalling (3.3) and (1.17) to treat the convex part β ε of F ε and the term involving N. At this point, by collecting everything, choosing δ small enough, and applying the Gronwall lemma, we obtain that In particular, this proves (3.36), so that T k = T .Then, by recalling (1.17)-(1.18)and the Poincaré inequality once more and rearranging, we conclude that Second a priori estimate.We now aim at recovering an estimate for the pressure p k .Thus, we construct for a.a.t ∈ (0, T ) and some constant C > 0. To this end, it suffices to apply Lemma 2.3 with an obvious choice of f and a.Then, we test (3.23),written at the time t, by q k (t).However, we avoid writing the time t for brevity.Recalling (1.6) and (3.24), we obtain that and the definition of q k , as well as the Hölder, Sobolev, and Young inequalities, yield V .Now, we rearrange, take the power of exponent 2/3, and integrate over (0, T ).By accounting for (3.44) and the Young inequality once more, we deduce that and it remains to estimate the last integral.By interpolation, we have the continuous embedding L ∞ (0, T ; H) ∩ L 2 (0, T ; L 6 (Ω)) ֒→ L 4 (0, T ; L 3 (Ω)) .
Since (3.44) implies that N σ (ϕ k , σ k ) is bounded in L ∞ (0, T ; H) ∩ L 2 (0, T ; V ) and the continuous embedding V ֒→ L 6 (Ω) holds, the integral at hand is uniformly bounded.Therefore, we have proved that By the way, the argument used for N σ (ϕ k , σ k ) also applies to σ k , so that Third a priori estimate.We test (3.26) by the admissible function −∆ϕ k and integrate in time.We obtain that As for the first term on the right-hand side, we recall (2.11) and Remark 3.3, to see that it is bounded.Since f is bounded in L 2 (0, T ; H) by (3.44), we deduce that the same holds for ∆ϕ k .Then, elliptic regularity yields that Fourth a priori estimate.We take any ζ ∈ L 4 (0, T ; V ) and define ζ k as follows: is the H-projection of ζ(t) onto V k for a.a.t ∈ (0, T ).Then, for a.a.t ∈ (0, T ), we test (3.27),written at the time t, by ζ k (t).However, we do not write the time t for simplicity.By also accounting for (3.24), we have that By the same remark, we have that . Hence, we infer that

Colli -Gilardi -Signori -Sprekels
To deal with the next term, we integrate by parts and obtain We can replace the L 6 norm by the V norm since V ֒→ L 6 (Ω), and the L 2 norm of ∇ζ k by ζ V .Moreover, since the two-dimensional embedding H 1/2 (Γ) ֒→ L 4 (Γ) holds true and the trace operator is continuous from V to H 1/2 (Γ), we can estimate the last product as follows: Similarly, we have that where, for clarity, we point out that σ Σ here means the norm of σ Σ (t) in L 2 (Γ).At this point, we collect all these inequalities and estimates and integrate over (0, T ).Omitting the integration variable t for brevity, we have that Therefore, by using the Hölder inequality, we have that Finally, we account for (3.44), (3.46), and Lemma 2.5, to conclude that This means that ∂ t σ k is bounded in the dual space of L 4 (0, T ; V ), i.e., that Fifth a priori estimate.We aim at estimating the time derivative ∂ t (εµ k + ϕ k ).To this end, we take any φ ∈ L 2 (0, T ; V 0 ) and define φ k : Q → R by setting for a.a.t ∈ (0, T ): is the H-projection of φ(t) onto V 0 k .Then, we rearrange (3.25), written at the time t, test it by φ k (t), and integrate over (0, T ).We obtain that For the first two terms on the right-hand side, we invoke the Lipschitz continuity of S ϕ and (3.44) to immediately obtain that while the last one needs some care.By the Hölder and Sobolev inequalities, (3.5) and (3.44) once more, we obtain that Since φ k L 2 (0,T ;V 0 ) ≤ φ L 2 (0,T ;V 0 ) by Remark 3.3 and we can replace φ by −φ, we infer that Unfortunately, just ∂ t µ k is V 0 k -valued while ∂ t ϕ k is not, so that we only have that However, on account of (3.44) and Remark 3.3 once more, we can write Combining this with (3.51) yields Passage to the limit as k → ∞.We collect the basic estimates we have proved, namely, (3.44), (3.45), (3.47), (3.50), and (3.52), and apply well-known weak and weak star compactness results.Since ε is fixed, for some (not relabeled) subsequence and suitable limit functions, we have, as k → ∞, p k → p ε weakly in L 4/3 (0, T ; H) , (3.54) In view of (3.56), note that (3.58) is equivalent to µ k → µ ε weakly in H 1 (0, T ; V * 0 ).Then, the quintuple (v ε , p ε , µ ε , ϕ ε , σ ε ) satisfies (3.7)- (3.11), as well as the initial conditions (3.17) (recall the Remarks 3.1 and 3.3), and we aim at proving that it solves the whole approximating problem.By linearity (cf.(1.18)), we deduce that, as k → ∞, On the other hand, the above convergence properties and the Lipschitz continuity imply that, as k → ∞, Now, we claim that Assume ζ k and ζ as said.From (3.57) and the strong compactness results already quoted, we deduce that, as k → ∞, From this, we deduce that, as k → ∞, owing to the interpolation theory in Hilbert spaces (see, e.g., [26,Ch. I]).We know that ∇ is continuous from V into H and from H into H −1 (Ω) = (H 1 0 (Ω)) * .Then, it is continuous from H 3/4 (Ω) = (V, H) 1/4 into (H, (H 1 0 (Ω)) * ) 1/4 .Since 1/4 < 1/2, the last space is (H Then, (3.66) follows.On the other hand, we have that where the last (three-dimensional) embedding follows from a particular case of the embedding Therefore, (3.65) holds true since L 2 (0, T ; H −1/2 (Ω)) = (L 2 (0, T ; H 1/2 (Ω))) * .By the same argument (by replacing the velocities by the N σ terms, essentially), one obtains that On account of the strong convergence in (3.69), recalling (3.58) and (3.62), and letting k tend to infinity, we conclude that By accounting for (3.65) and (3.69), we conclude that, as k → ∞, for all V -valued step functions z and ζ.Thus, (3.15)-(3.16)hold as well, and the proof of Theorem 3.2 is complete.

Existence of a Weak Solution
In this section, we prove Theorem 2.1.We start from the approximating problem analyzed in the previous section and let ε tend to zero.Since we did not prove uniqueness for the approximating solution, we take a particular one, namely, the solution we have constructed above.This ensures a number of bounds.Indeed, by the estimates established for the discrete solution and the semicontinuity of the norms, it is clear that for a positive constant c independent of ε.However, we need some additional estimates.
.40) At the same time, we test (3.25) and (3.26) by µ k and ∂ t ϕ k , respectively, and integrate with respect to time.We obtain that for some of the equations, and for simplicity we choose step functions for all of them).So, we assume that φ is a V 0 -valued step function and that z and ζ are V -valued step functions.However, we have to be careful, since we are starting from the discrete setting.andζhave a finite number of values.Next, we test(3.25),written at the time t, by φ k (t) and integrate over (0, T ).Upon rearranging, we obtain that So, given φ, z and ζ as said, we introduce φ k as we did in proving our fifth estimate, i.e., by defining φ k (t) as the H-projection of φ(t) onto V 0 k , for a.a.t ∈ (0, T ).Similarly, we define z k and ζ k starting from z and ζ, now with V k instead of V 0 k .By accounting for Remark 3.3, we point out at once that, as k → ∞,φ k → φ, z k → z,and ζ k → ζ, strongly in L ∞ (0, T ; V ), (3.69) since φ, z .70) What we have obtained is an integrated version of (3.14), which is equivalent to (3.14) itself since it holds for every V 0 -valued step function φ.Similarly, we test (3.26) and (3.27) written at the time t by z k (t) and ζ k (t), respectively, and integrate over (0, T ).