Global weak solutions and asymptotic limits of a Cahn--Hilliard--Darcy system modelling tumour growth

We study the existence of weak solutions to a Cahn--Hilliard--Darcy system coupled with a convection-reaction-diffusion equation through the fluxes, through the source terms and in Darcy's law. The system of equations arises from a mixture model for tumour growth accounting for transport mechanisms such as chemotaxis and active transport. We prove, via a Galerkin approximation, the existence of global weak solutions in two and three dimensions, along with new regularity results for the velocity field and for the pressure. Due to the coupling with the Darcy system, the time derivatives have lower regularity compared to systems without Darcy flow, but in the two dimensional case we employ a new regularity result for the velocity to obtain better integrability and temporal regularity for the time derivatives. Then, we deduce the global existence of weak solutions for two variants of the model; one where the velocity is zero and another where the chemotaxis and active transport mechanisms are absent.


Introduction
In recent years there has been an increased focus on the mathematical modelling and analysis of tumour growth. Many new models have been proposed and numerical simulations have been carried out to provide new and important insights on cancer research, see for instance [8] and [13,Chap. 3]. In this work we analyse a diffuse interface model proposed in [20], which models a mixture of tumour cells and healthy cells in the presence of an unspecified chemical species acting as a nutrient. More precisely, for a bounded domain Ω ⊂ R d where the cells reside and T > 0, we consider the following set of equations, in Ω × (0, T ) =: Q,

1d)
∂ t σ + div (σv) = div (n(ϕ)(D∇σ − χ∇ϕ)) − S in Q. (1.1e) Here, v denotes the volume-averaged velocity of the mixture, p denotes the pressure, σ denotes the concentration of the nutrient, ϕ ∈ [−1, 1] denotes the difference in volume fractions, with {ϕ = 1} representing the unmixed tumour tissue, and {ϕ = −1} representing the surrounding healthy tissue, and µ denotes the chemical potential for ϕ. The model treats the tumour and healthy cells as inertia-less fluids, leading to the appearance of a Darcy-type subsystem with a source term Γ v . The order parameter ϕ satisfies a convective Cahn-Hilliard type equation with additional source term Γ ϕ , and similarly, the nutrient concentration σ satisfies a convection-reaction-diffusion equation with a non-standard flux and a source term S. We refer the reader to [20, §2] for the derivation from thermodynamic principles, and to [20, §2.5] for a discussion regarding the choices for the source terms Γ ϕ , Γ v and S.
The positive constants K and D denote the permeability of the mixture and the diffusivity of the nutrient, m(ϕ) and n(ϕ) are positive mobilities for ϕ and σ, respectively. The parameter χ ≥ 0 regulates the chemotaxis effect (see [20] for more details), Ψ(·) is a potential with two equal minima at ±1, A and B denote two positive constants related to the thickness of the diffuse interface and the surface tension.
We supplement the above with the following boundary and initial conditions Here ϕ 0 , σ 0 and σ ∞ are given functions and b > 0 is a constant. We denote ∂ n f := ∇f · n as the normal derivative of f at the boundary ∂Ω, where n is the outer unit normal. Associated to (1.1) is the free energy density N (ϕ, σ) for the nutrient, which is defined as ∂ t ϕ + div (ϕv) = div (m(ϕ)∇µ) + Γ ϕ , (1.4c)

4d)
∂ t σ + div (σv) = div (n(ϕ)∇N ,σ ) − S, (1.4e) which is the general phase field model proposed in [20]. In this work we do not aim to analyse such a model with a general free energy density N (ϕ, σ), but we will focus solely on the choice (1.3) and the corresponding model (1.1)-(1.2). Our goal in this work is to prove the existence of weak solutions (see Definition 2.1 below) of (1.1)-(1.2) in two and three dimensions. Moreover, one might expect that by setting Γ v = 0 and then sending b → 0 and K → 0, the weak solutions to (1.1)-(1.2) will converge (in some appropriate sense) to the weak solutions of ∂ t ϕ = div (m(ϕ)∇µ) + Γ ϕ in Q, (1.5a) ∂ t σ = div (n(ϕ)(D∇σ − χ∇ϕ)) − S in Q, (1.5c) 0 = ∂ n ϕ = ∂ n µ = ∂ n σ on Σ. (1.5d) We denote (1.5) as the limit system of vanishing permeability, where the effects of the volume-averaged velocity are neglected. By substituting for some non-negative function f (ϕ) leads to the model derived in [21]. The specific choices for Γ ϕ and S in (1.6) are motivated by linear phenomenological laws for chemical reactions. The analysis of (1.5) with the parameters D = 1, χ = 0, n(ϕ) = m(ϕ) = 1 has been the subject of study in [5,7,6,16], where well-posedness and long-time behaviour have been established for a large class of functions Ψ(ϕ) and f (ϕ). Alternatively, one may consider the following choice of source terms where λ p , λ a , λ c are non-negative constants representing the tumour proliferation rate, the apoptosis rate, and the nutrient consumption rate, respectively, and h(ϕ) is a non-negative interpolation function such that h(−1) = 0 and h(1) = 1. The above choices for Γ ϕ and S are motivated from the modelling of processes experienced by a young tumour.
The well-posedness of model (1.5) with the choice (1.7) has been studied by the authors in [17] and [18] with the boundary conditions (1.2) (neglecting (1.2b)) in the former and for non-zero Dirichlet boundary conditions in the latter. It has been noted in [17] that the wellposedness result with the boundary conditions (1.2) requires Ψ to have at most quadratic growth, which is attributed to the presence of the source term Γ ϕ µ = h(ϕ)µ(λ p σ−λ a ) when deriving useful a priori estimates. Meanwhile in [18] the Dirichlet boundary conditions and the application of the Poincaré inequality allows us to overcome this restriction and allow for Ψ to be a regular potential with polynomial growth of order less than 6, and by a Yosida approximation, the case where Ψ is a singular potential is also covered.
We also mention the work of [19] that utilises a Schauder's fixed point argument to show existence of weak solutions for Ψ with quartic growth and Γ ϕ , S as in (1.7). This is based on first deducing that σ is bounded by a comparison principle, leading to Γ ϕ ∈ L ∞ (Ω). Then, the standard a priori estimates are derived for a Cahn-Hilliard equation with bounded source terms. The difference between [19] and [17,18] is the absence of the chemotaxis and active transport mechanisms, i.e., χ = 0, so that the comparison principle can be applied to the nutrient equation. We refer to [9] for the application of a similar procedure to a multi-species tumour model with logarithmic potentials.
On the other hand, by sending b → 0 and χ → 0 in (1.1), we should obtain weak ∂ t σ + div (σv) = div (n(ϕ)D∇σ) − S in Q, (1.8e) 0 = ∂ n ϕ = ∂ n µ = ∂ n σ = v · n on Σ. (1.8f) We denote (1.8) as the limit system of vanishing chemotaxis. If the source terms Γ v and Γ ϕ are independent of σ, then (1.8) consists of an independent Cahn-Hilliard-Darcy system and an equation for σ which is advected by the volume-averaged velocity field v. In the case where there is no nutrient and source terms, i.e., σ = Γ v = Γ ϕ = 0, global existence of weak solutions in two and three dimensions has been established in [14] via the convergence of a fully discrete and energy stable implicit finite element scheme. For the well-posedness and long-time behaviour of strong solutions, we refer to [25]. Meanwhile, in the case where Γ v = Γ ϕ is prescribed, global weak existence and local strong well-posedness for (1.8) without nutrient is shown in [22]. We also mention the work of [3] on the well-posedness and long-time behaviour of a related system also used in tumour growth, known as the Cahn-Hilliard-Brinkman system, where in (1.8) without nutrient an additional viscosity term is added to the left-hand side of the velocity equation (1.8b) and the mass exchange terms Γ v and Γ ϕ are set to zero. The well-posedness of a nonlocal variant of the Cahn-Hilliard-Brinkman system has been investigated in [10]. Furthermore, when K is a function depending on ϕ, the model (1.8) with σ = Γ v = Γ ϕ = 0 is also referred to as the Hele-Shaw-Cahn-Hilliard model (see [23,24]). In this setting, K(ϕ) represents the reciprocal of the viscosity of the fluid mixture. We refer to [30] concerning the strong well-posedness globally in time for two dimensions and locally in time for three dimensions when Ω is the d-dimensional torus. Global well-posedness in three dimensions under additional assumptions and long-time behaviour of solutions to the Hele-Shaw-Cahn-Hilliard model are investigated in [29].
We point out that from the derivation of (1.1) in [20], the source terms Γ v and Γ ϕ are connected in the sense that Γ v is related to sum of the mass exchange terms for the tumour and healthy cells, and Γ ϕ is related to the difference between the mass exchange terms. Thus, if Γ ϕ would depend on the primary variables ϕ, σ or µ, then one expects that Γ v will also depend on the primary variables. Here, we are able to prove existence of weak solutions for Γ ϕ of the form (2.1), which generalises the choices (1.6) and (1.7), but in exchange Γ v has to be considered as a prescribed function. This is attributed to the presence of the source term Γ v ϕµ + D 2 |σ| 2 when deriving useful a priori estimates. We see that if Γ v depends on the primary variables, we obtain triplet products which cannot be controlled by the usual regularity of ϕ, µ and σ in the absence of a priori estimates.
In this work we attempt to generalise the weak existence results for the models studied in [5,16,17,18,22,25] by proving that the weak solutions of (1.1) with Γ v = 0 converge (in some appropriate sense) to the weak solutions of (1.5) as b → 0 and K → 0, and the weak solutions of (1.1) converge to the weak solutions of (1.8) as b → 0 and χ → 0. This paper is organised as follows. In Section 2 we state the main assumptions and the main results. In Section 3 we introduce a Galerkin procedure and derive some a priori estimates for the Galerkin ansatz in Section 4 for the case of three dimensions. We then pass to the limit in Section 5 to deduce the existence result for three dimensions, while in Section 6 we investigate the asymptotic behaviour of solutions to (1.1) as K → 0 and χ → 0. In Section 7, we outline the a priori estimates for two dimensions and show that the weak solutions for two dimensions yields better temporal regularity than the weak solutions for three dimensions. In Section 8 we discuss some of the issues present in the analysis of (1.1) using different formulations of Darcy's law and the pressure, and with different boundary conditions for the velocity and the pressure.
Notation. For convenience, we will often use the notation L p := L p (Ω) and W k,p := W k,p (Ω) for any p ∈ [1, ∞], k > 0 to denote the standard Lebesgue spaces and Sobolev spaces equipped with the norms · L p and · W k,p . In the case p = 2 we use H k := W k,2 and the norm · H k . For the norms of Bochner spaces, we will use the notation L p (X) := L p (0, T ; X) for Banach space X and p ∈ [1, ∞]. Moreover, the dual space of a Banach space X will be denoted by X * , and the duality pairing between X and X * is denoted by ·, · X,X * . For d = 2 or 3, let H d−1 denote the (d − 1) dimensional Hausdorff measure on ∂Ω, and we denote R d -valued functions and any function spaces consisting of vector-valued/tensor-valued functions in boldface. We will use the notation Df to denote the weak derivative of the vector function f .
Useful preliminaries. For convenience, we recall the Poincaré inequality: There exists a positive constant C p depending only on Ω such that, for all f ∈ H 1 , where f := 1 |Ω| Ω f dx denotes the mean of f . The Gagliardo-Nirenberg interpolation inequality in dimension d is also useful (see [15,Thm. 10 Let Ω be a bounded domain with Lipschitz boundary, and f ∈ W m,r ∩ L q , 1 ≤ q, r ≤ ∞. For any integer j, 0 ≤ j < m, suppose there is α ∈ R such that Then, there exists a positive constant C depending only on Ω, m, j, q, r, and α such that (1.10) We will also use the following Gronwall inequality in integral form (see [17,Lem. 3.1] for a proof): Let α, β, u and v be real-valued functions defined on [0, T ]. Assume that α is integrable, β is non-negative and continuous, u is continuous, v is non-negative and integrable. If u and v satisfy the integral inequality then it holds that To analyse the Darcy system, we introduce the spaces Then, the Neumann-Laplacian operator −∆ N : H 1 ∩ L 2 0 → (H 1 ) * 0 is positively defined and self-adjoint. In particular, by the Lax-Milgram theorem and the Poincaré inequality (1.9) with zero mean, the inverse operator (−∆ N ) −1 : (H 1 ) * 0 → H 1 ∩ L 2 0 is well-defined, and we set u : −∆u = f in Ω, ∂ n u = 0 on ∂Ω.

Main results
We make the following assumptions.
for positive constants m 0 , m 1 , n 0 and n 1 .
(A3) Γ ϕ and S are of the form where Θ ϕ , Θ S : R 2 → R are continuous bounded functions with Θ ϕ non-negative, and Λ ϕ , Λ S : R 2 → R are continuous with linear growth for some positive constant R 0 .
(A5) Ψ ∈ C 2 (R) is a non-negative function satisfying and either one of the following,

1.
if Θ ϕ is non-negative and bounded, then if Θ ϕ is positive and bounded, that is,

7)
for some positive constants R 1 , R 2 , R 3 , R 4 , R 5 , R 6 . Furthermore we assume that The initial and boundary data satisfy We point out that some of the above assumptions are based on previous works on the well-posedness of Cahn-Hilliard systems for tumour growth. For instance, (2.5) and (2.8) reflect the situation encountered in [17], where if Θ ϕ = 0, i.e., Γ ϕ is independent of µ, then the derivation of the a priori estimate requires a quadratic potential. But in the case where (2.6) is satisfied, we can allow Ψ to be a regular potential with polynomial growth of order less than 6, and by a Yosida approximation, we can extend our existence results to the situation where Ψ is a singular potential, see for instance [18]. Moreover, the condition (2.8) is a technical assumption based on the fact that the second term of the nutrient free energy χσ(1 − ϕ) does not have a positive sign.
Meanwhile, the linearity of the source terms Γ ϕ and S with respect to the chemical potential µ assumed in (2.1) is a technical assumption based on the expectation that, at best, we have weak convergence for Galerkin approximation to µ, which is in contrast with ϕ and σ where we might expect a.e convergence and strong convergence for the Galerkin approximations. Moreover, if we consider for a non-negative function f (ϕ), then we obtain the source terms in [5,16,21].
Compared to the set-up in [22], in (A4) we prescribe a higher temporal regularity for the prescribed source term Γ v . This is needed when we estimate the source term Γ v D 2 |σ| 2 in the absence of a priori estimates, see Section 4.1.2 for more details. The mean zero condition is a consequence of the no-flux boundary condition v · n = 0 on ∂Ω and the divergence equation (1.1a). In particular, we can express the Darcy subsystem (1.1a)-(1.1b) as an elliptic equation for the pressure p: Solutions to (2.9) are uniquely determined up to an arbitrary additive function that may only depend on time, and thus without loss of generality, we impose the condition p =

1
|Ω| Ω p dx = 0 to (2.9). We may then define p as Remark 2.1. In the case Γ v = 0, one can also consider the assumption instead of (2.6), which holds automatically if Γ ϕ and S are chosen to be of the form (1.6).
In fact this property is used in [5,16].
We make the following definition.
Neglecting the nutrient σ, we observe that our choice of function spaces for (ϕ, µ, p, v) coincide with those in [22,Defn. 2.1(i)]. In contrast to the usual L 2 (0, T ; (H 1 ) * )-regularity (see [5,16]) we obtain a less regular time derivative ∂ t ϕ. The drop in the time regularity from 2 to 8 5 is attributed to the convection term div (ϕv) belonging to L 8 5 (0, T ; (H 1 ) * ). The same is true for the regularity for the time derivative ∂ t σ in L 5 4 (0, T ; (W 1,5 ) * ) as the convection term div (σv) lies in the same space. We refer the reader to the end of Section 4.3 for a calculation motivating the choice of function spaces for div (σv) and ∂ t σ. Furthermore, the embedding of L ∞ (0, T ; [28,§8,Cor. 4] guarantees that the initial condition for ϕ is meaningful. However, for σ we have the embedding L ∞ (0, T ; L 2 ) ∩ W 1, 5 4 (0, T ; (W 1,5 ) * ) ⊂⊂ C 0 ([0, T ]; (H 1 ) * ), and so σ(0) makes sense as a function in (H 1 ) * . Thus, the initial condition σ 0 is attained as an equality in (H 1 ) * . We now state the existence result for (1.1)-(1.2). and in addition satisfies (2.14) where the constant C does not depend on (ϕ, µ, σ, v, p) and is uniformly bounded for b, χ ∈ (0, 1] and is also uniformly bounded for K ∈ (0, 1] when Γ v = 0. The regularity result (2.13) is new compared to estimates for weak solutions in [22], which arises from a deeper study of the Darcy subsystem, and can be obtained even in the absence of the nutrient. We mention that higher regularity estimates for the pressure p in L 2 (0, T ; H 2 ) and the velocity v in L 2 (0, T ; H 1 ) are also established in [22], but these are for strong solutions local in time in three dimensions and global in time for two dimensions.
We now investigate the situation in two dimensions, where the Sobolev embeddings in two dimensions yields better integrability exponents.
The proof of Theorem 2.2 is similar to that of Theorem 2.1, and hence the details are omitted. In Section 7 we will only present the derivation of a priori estimates. It is due to the better exponents for embeddings in two dimensions and the regularity result for the velocity that we obtain better regularities for the time derivatives ∂ t ϕ and ∂ t σ, namely ∂ t σ(t) belongs to the dual space (H 1 ) * for a.e. t ∈ (0, T ). Furthermore, as mentioned in Remark 7.1 below, if we only have v ∈ L 2 (0, T ; L 2 ), then the convection term div (σv) and the time derivative ∂ t σ would only belong to the dual space L 4 3 (0, T ; (W 1,4 ) * ). However, even with the improved temporal regularity, as ∂ t σ / ∈ L 2 (0, T ; (H 1 ) * ), we do not have a continuous embedding into the space C 0 ([0, T ]; L 2 ) and so σ(0) may not be well-defined as an element of L 2 .
We now state the two asymptotic limits of (1.1) for three dimensions, and note that analogous asymptotic limits also hold for two dimensions.

Galerkin approximation
We will employ a Galerkin approximation similar to the one used in [22]. For the approximation, we use the eigenfunctions of the Neumann-Laplacian operator {w i } i∈N . Recall that the inverse Neumann-Laplacian operator L : Furthermore, let {f n } n∈N ⊂ L 2 0 denote a sequence with corresponding solution sequence {z n = Lf n } n∈N ⊂ H 1 ∩ L 2 0 . By elliptic regularity theory, we have that z n ∈ H 2 N for all n ∈ N. Then, by reflexive weak compactness theorem and Rellich-Kondrachov theorem, there exists a subsequence such that z n j → z ∈ H 1 ∩ L 2 0 as j → ∞. Thus, by the spectral theorem, the operator L admits a countable set of eigenfunctions {v n } n∈N that forms a complete orthonormal system in L 2 0 . The eigenfunctions of the Neumann-Laplacian operator is then given by Elliptic regularity theory gives that w i ∈ H 2 N and for every g ∈ H 2 N , we obtain for where λ i is the corresponding eigenvalue to w i . This shows that ∆g k converges strongly to ∆g in L 2 . Making use of elliptic regularity theory again gives that g k converges strongly to g in H 2 N . Thus the eigenfunction {w i } i∈N of the Neumann-Laplace operator forms an orthonormal basis of L 2 and is also a basis of H 2 N . Later in Section 5, we will need to use the property that H 2 N is dense in H 1 and W 1,5 . We now sketch the argument for the denseness of H 2 N in W 1,5 and the argument for H 1 follows in a similar fashion.
Proof. Take g ∈ W 1,5 , as Ω has a C 3 -boundary, by standard results [12, Thm. 3, §5.3.3] there exists a sequence g n ∈ C ∞ (Ω) such that g n → g strongly in W 1,5 . Let ε > 0 be fixed, and define D ε := {x ∈ Ω : dist(x, ∂Ω) ≤ ε}. Let ζ ε ∈ C ∞ c (Ω) be a smooth cut-off function such that ζ ε = 1 in Ω \ D ε and ζ ε = 0 in D ε 2 . As g n ∈ C ∞ (Ω), its trace on ∂Ω is well-defined. Choosing ε sufficiently small allows us to use a classical result from differential geometry about tubular neighbourhoods, i.e., for any z ∈ Tub ε (∂Ω) := {x ∈ R d : |dist(z, ∂Ω)| ≤ ε} there exists a unique y ∈ ∂Ω such that where n is the outer unit normal of ∂Ω. We consider a bounded smooth function f n,ε : We now define the smooth function G n,ε as By construction, the values of the function f n,ε in D ε ⊂ Tub ε (∂Ω) are constant in the normal direction, so ∇G n,ε · n = 0 on ∂Ω and thus G n,ε ∈ H 2 N . Furthermore, we compute that Using that g n , f n,ε are smooth functions on Ω and that the Lebesgue measure of D ε tends to zero as ε → 0 we have the strong convergence of G n,ε to g n in L 5 . For the difference in the gradients, we use that ζ ε → 1 a.e. in Ω, Lebesgue's dominated convergence theorem and the boundedness of ∇g n and ∇f n,ε to deduce that For the remaining term (g n − f n,ε )∇ζ ε L 5 we use that the support of ∇ζ ε lies in D ε \ D ε 2 and for any z ∈ D ε \ D ε 2 , |f n,ε (z) − g n (z)| = |g n (y) − g n (y + dist(z, ∂Ω)n(y))| That is, f n,ε converges uniformly to This shows that G n,ε converges strongly to g n in W 1,5 .
We denote W k := span{w 1 , . . . , w k } as the finite dimensional space spanned by the first k basis functions and consider and the following Galerkin ansatz: For 1 ≤ j ≤ k, where we define the Galerkin ansatz for the pressure p k and the velocity field v k by and we set Note that in (3.3), the properties Γ v ∈ L 2 0 and ∇ϕ k · n = 0 on ∂Ω show that the term inside the bracket belongs to L 2 0 and hence p k is well-defined. Let M and S denote the following mass and stiffness matrices, respectively: Thanks to the orthonormality of {w i } i∈N in L 2 , we see that M is the identity matrix. It is convenient to define the following matrices with components as the corresponding vectors, so that we obtain an initial value problem for a system of equations for α k : and we supplement (3.5) with the initial conditions for some constant C not depending on k.
We can substitute (3.5b), (3.5d) and (3.5e) into (3.5a) and (3.5c), and obtain a coupled system of ordinary differential equations for α k and γ k , where S k m , C k and S k n depend on the solutions α k and γ k in a non-linear manner. Continuity of m(·), n(·), Ψ ′ (·) and the source terms, and the stability of (−∆ N ) −1 under perturbations imply that the right-hand sides of (3.5) depend continuously on (α k , γ k ). Thus, we can appeal to the theory of ODEs (via the Cauchy-Peano theorem [4, Chap. 1, Thm. 1.2]) to infer that the initial value problem (3.5)-(3.6) has at least one local solution pair (α k , γ k ) defined on [0, t k ] for each k ∈ N.
We may define β k via the relation (3.5b) and hence the Galerkin ansatz ϕ k , µ k and σ k can be constructed from (3.1). Then, we can define p k and v k via (3.3) and (3.4), respectively. Furthermore, as the basis function w j belongs to H 2 for each j ∈ N, by the Sobolev embedding H 2 ⊂ L ∞ , we obtain that div (w i ∇w j ) ∈ L 2 for i, j ∈ N and hence the function div ((µ k + χσ k )∇ϕ k ) belongs to L 2 . Then, by elliptic regularity theory, we Next, we show that the Galerkin ansatz can be extended to the interval [0, T ] using a priori estimates.

A priori estimates
In this section, the positive constants C are independent of k, Γ v , K, b and χ, and may change from line to line. We will denote positive constants that are uniformly bounded for b, χ ∈ (0, 1] and are also uniformly bounded for K ∈ (0, 1] when Γ v = 0 by the symbol E.
We first state the energy identity satisfied by the Galerkin ansatz. Let δ ij denote the Kronecker delta.
, and then summing the product from j = 1 to k lead to Here, we used that w 1 = 1 and ∇w 1 = 0. Then, summing the three equations leads to Next, multiplying (3.4) with 1 K v k , integrating over Ω and integrating by parts gives where we used that div v k = Γ v and v k · n = 0 on ∂Ω. Similarly, we see that In particular, we have To derive the first a priori estimate for the Galerkin ansatz, it suffices to bring (4.4) into a form where we can apply Gronwall's inequality. We start with estimating the boundary term on the right-hand side of (4.4). By Hölder's inequality and Young's inequality, By the trace theorem and the growth condition (2.4), we have where the positive constant C tr from the trace theorem only depends on Ω, and so (4.6)

Estimation of the source terms
For the source term that appears on the right-hand side of (4.4) we will divide its analysis into two parts. We first analyse the part involving Γ v , which will involve a closer look at the Darcy subsystem to deduce an estimate on p k L 2 . For the remainder Γ ϕ,k µ k − S k N k ,σ term we will estimate it differently based on the assumptions on Θ ϕ .

Pressure estimates
Before we estimate the source terms involving Γ v , we look at the Darcy subsystem, which can be expressed as an elliptic equation for the pressure (we will drop the subscript k for clarity) The following lemma is similar to [22,Lem. 3.1], and the hypothesis is fulfilled by the Galerkin ansatz.
Let Ω ⊂ R 3 be a bounded domain with C 3 -boundary. Given ϕ ∈ H 2 N , µ, σ ∈ H 1 , the source term Γ v ∈ L 2 0 , and the function p satisfying the above elliptic equation (4.7). Then, the following estimate hold Then, testing with f and integrating over Ω, applying integration by parts and the Poincaré inequality (1.9) leads to for positive constants c and C depending only on C p . Elliptic regularity theory then gives with a positive constant C depending only on Ω. Returning to the pressure system, we observe from (2.10) and the above that for some positive constant C depending only on C p . Note that the third term on the right-hand side can be estimated as We now consider estimating the second term on the right-hand side of (4.12). By assumption µ, σ ∈ H 1 and ϕ ∈ H 2 N , we have that and so if we consider the function h := (−∆ N ) −1 ( div ((µ − µ + χσ)∇ϕ)), then we obtain that Ω ∇h · ∇ζ dx = Ω −(µ − µ + χσ)∇ϕ · ∇ζ dx ∀ζ ∈ H 1 (4.15) must hold, and by (4.14) and the Poincaré inequality (1.9) with zero mean it holds that h ∈ H 1 ∩ L 2 0 . We now define f := (−∆ N ) −1 (h) ∈ H 2 N , and consider testing with ζ = f in (4.15), leading to Since f ∈ H 2 N , elliptic regularity theory and Hölder's inequality gives where the constant C depends on Ω and the constant in (4.11). Thus we obtain  for some constant C depending only on Ω. By the Sobolev embedding H 1 ⊂ L 6 (with constant C Sob that depends only on Ω) and the Poincaré inequality, we find that (4.17) Substituting the above elements into (4.12) yields (4.8).
Remark 4.1. We choose not to use the estimate obtained from substituting ζ = h in (4.15), where c is a positive constant depending only on C p , since by (4.14) we require control of ∇ϕ in the L 3 (Ω)-norm and this is not available when deriving the first a priori estimate. Thus, we make use of the auxiliary problem f = (−∆ N ) −1 (h) to derive another estimate on h L 2 that involves controlling ∇ϕ in the weaker L 3 2 (Ω)-norm.
Next, we state regularity estimates for the pressure and the velocity field. The hypothesis will be fulfilled for the Galerkin ansatz once we derived the a priori estimates in Section 4. Note that in Lemma 4.2 below, we consider a source term Γ v ∈ L 2 (0, T ; L 2 0 ), so that our new regularity results for the pressure and the velocity is also applicable to the setting considered in [22].
, and the function p satisfying (4.7). Then, for some positive constant C 1 depending only on Ω, and for some positive constant C 3 depending only on Ω.
Proof. From (4.7) we see that p satisfies p = 0 and Testing with ζ = p and applying the Hölder's inequality and the Poincaré inequality (1.9) gives Applying Hölder's inequality and the Sobolev embedding H 1 ⊂ L 6 yields that By the Gagliardo-Nirenberg inequality (1.10) with parameters j = 0, p = 3, r = 2, m = 2, d = 3 and q = 2, where C > 0 is a constant depending only on Ω. Then, the boundedness of µ, σ in L 2 (0, T ; H 1 ) and ϕ in L 2 (0, T ; By (4.22) we find that where the positive constant C depends only on Ω. As p = 0, by the Poincaré inequality (1.9), we see that for some positive constant C depending only on Ω. Next, we see that By the Gagliardo-Nirenberg inequality (1.10), we find that  and so, we have div ((µ + χσ)∇ϕ) L 2 ≤ C µ + χσ H 1 ϕ (4.26) That is, div ((µ + χσ)∇ϕ) ∈ L 2 . Since by assumption Γ v ∈ L 2 0 , using elliptic regularity theory, we find that p(t, ·) ∈ H 2 for a.e. t and there exists a constant C depending only on Ω, such that (4.27) Furthermore, from (4.26), we see that  For the velocity field v we estimate as follows. Let 1 ≤ i, j ≤ 3 be fixed, we obtain from (4.25), (4.29) Applying the same calculation as in (4.28) yields for some positive constant C depending only on Ω.

Source term from the Darcy system
To estimate the third source term of the energy equality we use Hölder's inequality to obtain By the Gagliardo-Nirenberg inequality (1.10) with j = 0, r = 2, m = 1, p = 4, q = 2 and α = 3 4 , we have By Young's inequality with Hölder exponents (i.e., ab ≤ ε p a p + ε −q/p q b q for 1 p + 1 q = 1 and ε > 0), we find that for some positive constant C depending only in n 0 , D and Ω. Then, by (4.17) we have where the positive constant C depends only on Ω, m 0 , n 0 and D. Here we point out that the assumption Γ v ∈ L 4 (0, T ; L 2 0 ) is needed. For the remainder term Γ v (p k − µ k ϕ k ), we find that where we used Then, by Applying the calculations in the proof of Lemma 4.1 (specifically (4.10), (4.16) and (4.17)), Hölder's inequality and Young's inequality, we find that where C is a positive constant depending only on |Ω|, C p , C Sob , D, n 0 and m 0 . Here we point out that if we applied (4.18) instead of (4.8) then we obtain a term containing ∇ϕ L 3 on the right-hand side and this cannot be controlled by the left-hand side of (4.4). Using (2.4) we have Then, we obtain the following estimate for some positive constant C depending only on R 1 , R 2 , Ω, m 0 , n 0 and D. Here we point out that it is crucial for the source term Γ v to be prescribed and is not a function of ϕ, µ and σ, otherwise the product term Γ v 4 L 2 σ k 2 L 2 and Γ v 2 L 2 ∇ϕ k 2 L 2 cannot be controlled in the absence of any a priori estimates. For the remaining source term Ω Γ ϕ,k µ k − S k N k ,σ dx we split the analysis into two cases and combine with (4.31) to derive an energy inequality.

Energy inequality for non-negative Θ ϕ
Suppose Θ ϕ is non-negative and bounded, and Ψ is a potential that satisfies (2.5). We will estimate the mean of µ k by setting j = 1 in (3.2b), and using the growth condition (2.5) to obtain Then, by the Poincaré inequality (1.9) and the growth condition (2.4), we find that (4.32) Note that by the specific form (2.1) for Γ ϕ we have that Moving the non-negative term Θ ϕ (ϕ k , σ k ) |µ k | 2 to the left-hand side of (4.4) and subsequently neglecting it, we estimate the remainder using the growth condition (2.3) and Hölder's inequality as follows (here we use the notation Λ ϕ,k := Λ ϕ (ϕ k , σ k )), where C is a positive constant depending only on R 0 and |Ω|. By Young's inequality, (4.32) and (4.30), we have 34) for some positive constant C depending only on |Ω|, R 0 , R 1 , R 2 , R 4 , A, D, C p and m 0 . Using the fact that we now estimate the right-hand side of (4.4) using (4.6), (4.31) and (4.34), which leads to for some positive constant C not depending on Γ v , K, b and χ. Integrating (4.35) with respect to t from 0 to s ∈ (0, T ] leads to for some positive constant C independent of Γ v , K, χ and b. Here we used σ 0 ∈ L 2 and ϕ 0 ∈ H 1 , which implies by the growth condition (2.5) that Ψ(ϕ 0 ) ∈ L 1 . Next, by Hölder's inequality and Young's inequality we have (4.37) Substituting (4.37) into (4.36) then yields for some positive constant C independent of Γ v , K, b and χ. Setting where we recall that E denotes a constant that is uniformly bounded for b, χ ∈ (0, 1] and is also uniformly bounded for K ∈ (0, 1] when Γ v = 0.

Energy inequality for positive Θ ϕ
Suppose Θ ϕ satisfies (2.6) and Ψ is a potential satisfying the growth condition (2.7). Similar to the previous case, we see that the specific form for Γ ϕ leads to We move the term Θ ϕ (ϕ k , σ k ) |µ k | 2 to the left-hand side of (4.4) and estimate the remainder as in (4.33). Using Young's inequality differently and also (4.30), we have , (4.41) for some positive constant C depending only on |Ω|, R 5 , R 1 , R 2 , A, D, and C p . Using (4.6), (4.31), (4.41) and the lower bound Θ ϕ ≥ R 5 , instead of (4.35) we obtain from (4.4) for some positive constant C independent of Γ v , K, b and χ. We point out the main difference between (4.35) and the above is the appearance of the term R 5 2 µ k 2 L 2 on the left-hand side. The positivity of Θ ϕ allows us to absorb the µ k 2 L 2 term on the right-hand side of (4.41) and thus we do not need to use (4.32), which was the main reason why Ψ has to be a quadratic potential for a non-negative Θ ϕ . Then, applying a similar argument as in Section 4.1.3, we arrive at an analogous energy inequality to (4.40), This a priori estimate implies that the Galerkin ansatz ϕ k , µ k , σ k and v k can be extended to the interval [0, T ]. To determine if p k can also be extended to the interval [0, T ] we require some higher order estimates for ϕ k in order to use (4.19).

Higher order estimates
Let Π k denote the orthogonal projection onto the finite-dimensional subspace W k . From (3.2b) we may view ϕ k as the solution to the following elliptic equation For the case where Ψ satisfies (2.5), as {ϕ k } k∈N is bounded in L ∞ (0, T ; H 1 ), we have that {Ψ ′ (ϕ k )} k∈N is also bounded in L ∞ (0, T ; H 1 ). Using the fact that our basis functions {w i } i∈N are the eigenfunctions of the inverse Neumann-Laplacian operator and is therefore orthogonal in H 1 , and the Sobolev embedding H 1 ⊂ L r for r ∈ [1,6], there exists a positive constant C independent of ϕ k such that Then, this implies that {Π k (Ψ ′ (ϕ k ))} k∈N is also bounded in L ∞ (0, T ; H 1 ). As the righthand side of (4.45a) belongs to H 1 for a.e. t ∈ (0, T ), and the boundary ∂Ω is C 3 , by elliptic regularity theory, we have For fixed m ∈ [1, 5), we define a sequence of positive numbers {l j } j∈N by It can be shown that {l j } j∈N is a strictly increasing sequence such that l j → ∞ as j → ∞. The Gagliardo-Nirenberg inequality (1.10) then yields the following continuous embedding At the first step, the boundedness of {ϕ k } k∈N in L ∞ (0, T ; H 1 ) yields which implies that {Π k (Ψ ′ (ϕ k ))} k∈N is bounded in L 2 (0, T ; L l 1 ). As the other terms on the right-hand side of (4.45) are bounded in L 2 (0, T ; H 1 ), elliptic regularity then yields that {ϕ k } k∈N is bounded in L 2 (0, T ; W 2,l 1 ), and thus in L 2m (0, T ; L ml 2 ) by (4.49). At the j-th step, we have {ϕ k } k∈N is bounded in L 2 (0, T, W 2,l j ) ∩ L 2m (0, T ; L ml j+1 ). Then, it holds that and so {Π k (Ψ ′ (ϕ k ))} k∈N is bounded in L 2 (0, T ; L l j+1 ). Elliptic regularity then implies that {ϕ k } k∈N is bounded in L 2 (0, T ; W 2,l j+1 ). We terminate the bootstrapping procedure once l j ≥ 6 for some j ∈ N. This occurs after a finite number of steps as lim j→∞ l j = ∞. Altogether, we obtain that {ϕ k } k∈N is bounded in L 2 (0, T ; W 2,6 ). From (4.48) it holds that and by the following continuous embeddings obtain from the Gagliardo-Nirenberg inequality (1.10), we find that {Π k (Ψ ′ (ϕ k ))} k∈N is bounded in L 2 (0, T ; H 1 ). Applying elliptic regularity once more leads to the boundedness of {ϕ k } k∈N in L 2 (0, T ; H 3 ). Consequently, the hypotheses of Lemma 4.2 are satisfied and we obtain that which implies that the Galerkin ansatz p k can be extended to the interval [0, T ].
Note that the temporal exponent decreases while the spatial exponent increases as r increases, and they intersect at the point r = 10 3 .
By the continuity of Λ ϕ we have The same arguments can be applied for the source term S and for the derivative Ψ ′ (ϕ) satisfying the linear growth condition (2.5). For potentials satisfying the growth condition (2.7), we refer to the argument in [18, §3.1.2].

Existence in two dimensions
We first derive an analogous result to Lemma 4.2 for two dimensions.

Discussion
Reformulations of Darcy's law and the pressure. Associated to Darcy's law (1.1b) is the term λ v := p − µϕ − D 2 |σ| 2 which will contribute the source term Γ v λ v in the energy identity (4.4). In [20, Rmk. 2.1] three other reformulations of Darcy's law (1.1b) and the pressure are considered: (R1) Let q := p − AΨ(ϕ) − B 2 |∇ϕ| 2 so that From the viewpoint of estimating the source term Γ v λ v , we see that (8.3a) has the advantage of being the simplest. Meanwhile, for (8.2a) the analysis for Γ v λ v is similar to that performed in Section 4.1.2, but for (8.1a) the main difficulty will be to estimate the terms (AΨ(ϕ) + B 2 |∇ϕ| 2 )Γ v and (− D 2 |σ| 2 + χσϕ)Γ v , which at first glance would require the assumption that Γ v ∈ L ∞ (Q), and obtaining an L 2 -estimate for the pressure q from the Darcy law (8.1b) would be difficult due to the term div (∇ϕ ⊗ ∇ϕ).
Other boundary conditions for the pressure and velocity. In [20, §2.4.4] the authors have discussed possible boundary conditions for the velocity and for the pressure. As discussed in Section 2 following Assumption 2.1, we require the source term Γ v to have zero mean due to the no-flux boundary condition v · n = 0 on ∂Ω. The general energy identity (with homogeneous Neumann boundary conditions for ϕ and µ) from [20, Equ. which motivates the consideration of a Robin-type boundary condition forp: g = ap − v · n = ap + K∂ np − K(Dσ + χ(1 − ϕ))∂ n σ on ∂Ω, for some given datum g and positive constant a. On one hand, this would allow us to consider source terms Γ v that need not have zero mean, but on the other hand, the analysis of the Darcy system becomes more complicated. In particular, the weak formulation of the pressure system now reads as Ω K∇p · ∇ζ dx + ∂Ω apζ dH d−1 = Ω Γ v ζ + K (µ∇ϕ + (Dσ + χ(1 − ϕ)) ∇σ) · ∇ζ dx + ∂Ω gζ dH d−1 , and we observe that the term Dσ∇σ on the right-hand side belongs to L 1 as σ has at most H 1 -spatial regularity from the energy identity. Thus, it is not clear if the pressure system can be solved with the regularities stated in Lemma 4.1. A deeper study into the theory of linear elliptic equations with right-hand sides of the form div f where f ∈ L 1 is required.