Finite‐volume scheme for a degenerate cross‐diffusion model motivated from ion transport

An implicit Euler finite‐volume scheme for a degenerate cross‐diffusion system describing the ion transport through biological membranes is proposed. The strongly coupled equations for the ion concentrations include drift terms involving the electric potential, which is coupled to the concentrations through the Poisson equation. The cross‐diffusion system possesses a formal gradient‐flow structure revealing nonstandard degeneracies, which lead to considerable mathematical difficulties. The finite‐volume scheme is based on two‐point flux approximations with “double” upwind mobilities. The existence of solutions to the fully discrete scheme is proved. When the particles are not distinguishable and the dynamics is driven by cross diffusion only, it is shown that the scheme preserves the structure of the equations like nonnegativity, upper bounds, and entropy dissipation. The degeneracy is overcome by proving a new discrete Aubin–Lions lemma of “degenerate” type. Numerical simulations of a calcium‐selective ion channel in two space dimensions show that the scheme is efficient even in the general case of ion transport.


INTRODUCTION
The ion transport through biological channels plays an important role in all living organisms. On a macroscopic level, the transport can be described by nonlinear partial differential equations for the ion concentrations (or, more precisely, volume fractions) and the surrounding electric potential. A classical model for ion transport are the Poisson-Nernst-Planck equations [1], which satisfy Fick's law for the fluxes. However, this approach does not include size exclusion effects in narrow ion channels. Taking into account the finite size of the ions, one can derive from an on-lattice model in the diffusion limit another set of differential equations with fluxes depending on the gradients of all species [2,3]. These nonlinear cross-diffusion terms are common in multicomponent systems [4,Chap. 4].
In the general case, the evolution of the concentrations u i and fluxes  i of the ith ion species is governed by the equations for i = 1, …, n, where u 0 = 1 − ∑ n i=1 u i is the concentration (volume fraction) of the electro-neutral solvent, D i > 0 is a diffusion coefficient, > 0 is the (scaled) inverse thermal voltage, and z i ∈ R the charge of the ith species. Observe that we assumed Einstein's relation which says that the quotient of the diffusion and mobility coefficients is constant, and we call this constant 1/ . The electric potential is determined by the Poisson equation where 2 is the (scaled) permittivity constant and f = f (x) is a permanent background charge density.
Equations (1) and (2) are solved in a bounded domain Ω ⊂ R d . In order to match experimental conditions, the boundary Ω is supposed to consist of two parts, the insulating part Γ N , on which no-flux boundary conditions are prescribed, and the union Γ D of boundary contacts with external reservoirs, on which the concentrations are fixed. The electric potential is prescribed at the electrodes on Γ D . This leads to the mixed Dirichlet-Neumann boundary conditions where the boundary data (u i ) 1≤i≤n and Φ can be defined on the whole domain Ω. Finally, we prescribe the initial conditions The main mathematical difficulties of Equations (1) are the strong coupling and the fact that the diffusion matrix (A ij (u)), defined by A ij (u) = D i u i for i ≠ j and A ii (u) = D i (u 0 + u i ) is not symmetric and not positive definite. It was shown in Burger and coworkers [2,5] that system (1) possesses a formal gradient-flow structure. This means that there exists a (relative) entropy functional H[u] = ∫ Ω h(u)dx with the entropy density where u = (u 1 , …, u n ) and u 0 = 1 − ∑ n i=1 u i , such that (1) can be formally written as where B ii = D i u 0 u i , B ij = 0 for i ≠ j provide a diagonal positive definite matrix, and w j are the entropy variables, defined by We refer to [6,Lem. 7] for the computation of h/ u i . The entropy structure of (1) is useful for two reasons. First, it leads to L ∞ bounds for the concentrations. Indeed, the transformation (u, Φ)  → w to entropy variables can be inverted, giving u = u(w, Φ) with Then u i is positive and bounded from above, that is, This yields L ∞ bounds without the use of a maximum principle. Second, the entropy structure leads to gradient estimates via the entropy inequality where the constant C > 0 depends on the Dirichlet boundary data. Because of we achieve gradient estimates for u 1∕2 0 u i and u 1∕2 0 . Since u 0 may vanish locally, this does not give gradient bounds for u i , which expresses the degenerate nature of the cross-diffusion system. As a consequence, the flux has to be formulated in the terms of gradients of u 1∕2 0 u i and u 1∕2 0 only, namely The challenge is to derive a discrete version of this formulation. It turns out that (24) below is the right formulation in our context (assuming vanishing drift parts).
Our aim is to design a numerical approximation of (1) which preserves the structural properties of the continuous equations under simplifying assumptions (the coefficients D i are the same and the drift term vanishes). This suggests to use the entropy variables as the unknowns, as it was done in our previous work [6] with simulations in one space dimension. Unfortunately, we have not been able to perform a numerical convergence analysis with these variables. The reason is that we need discrete chain rules in order to formulate (7) on the discrete level and these discrete chain rules seem to be not easily available. Therefore, we use the original variables u i for the numerical discretization.
We propose a backward Euler scheme in time and a finite-volume scheme in space, based on two-point approximations. The key observation for the numerical discretization is that the fluxes can be written on each cell in a "double" drift-diffusion form, that is, both where v is the diffusion term and vF is the drift term. We discretize  and V by using a two-point flux aproximation with "double" upwind mobilities.
Under certain assumptions, the structure of the equations is preserved on the discrete level. Because of the drift-diffusion structure, we are able to prove that the scheme preserves the nonnegativity, which follows from a discrete minimum principle argument. It is well known that the maximum priciple generally does not hold for systems of equations. Therefore, it is not a surprise that the upper bound comes only at a price: We need to assume that the all diffusion coefficients D i are the same. Under this assumption, u 0 = 1 − ∑ n i=1 u i solves a drift-diffusion equation for which the (discrete) maximum principle can be applied. In order to prove that the scheme satisfies an entropy-dissipation inequality and also to complete the convergence analysis of the scheme successfully, we need a stronger additional assumption: We assume that the drift terms, and therefore the coupling with the Poisson equation, can be neglected. This means that our main results are obtained for a simplified degenerate cross-diffusion system, no more corresponding to the initial ion transport model but still of mathematical interest. Nevertheless, the scheme we propose can be applied to the full ion transport model, and this is done in the last section of this paper.
Our analytical results are stated and proved for no-flux boundary conditions on Ω. Mixed Dirichlet-Neumann boundary conditions could be prescribed as well, but the proofs would become even more technical. The main results are as follows.
• If D i = D for all i, we prove the existence of solutions to the fully discrete numerical scheme (theorem 1). If additionally the drift part vanishes, the solution is unique. The existence proof uses a topological degree argument in finite space dimensions, while the uniqueness proof is based on the entropy method of Gajewski [7], recently extended to cross-diffusion systems [6,8]. • If D i = D for all i, the scheme preserves the nonnegativity and upper bound for the concentrations. If additionally the drift part vanishes, convexity arguments show that the discrete entropy is dissipated with a discrete entropy production analogous to (7) (theorem 2). The assumption on vanishing drift terms is needed, since a discrete version of the sum ∑ n i=1 u i has to be controlled from below; see the discussion after theorem 2.
• If D i = D for all i and the drift part vanishes, the discrete solution converges to a continuous solution to (1) as the mesh size tends to zero (theorem 3). The proof is based on a priori estimates obtained from the discrete entropy inequality. The compactness is derived from a new discrete Aubin-Lions lemma, which takes into account the nonstandard degeneracy of the equations; see lemma 10 in the appendix. • Numerical experiments for a calcium-selective ion channel in two space dimensions show the dynamical behavior of the solutions and their large-time asymptotics to the equilibrium. The tests indicate that the order of convergence in the L 1 norm is one.
In the literature, there exist some results on finite-volume schemes for cross-diffusion systems. An upwind two-point flux approximation similar to our discretization was recently used in Ait Hammou Oulhaj [9] for a seawater intrusion cross-diffusion model. A two-point flux approximation with a nonlinear positivity-preserving approximation of the cross-diffusion coefficients, modeling the segregation of a two-species population, was suggested in Andreianov and coworkers [10], assuming positive definiteness of the diffusion matrix. The Laplacian structure of the population model (still for positive definite matrices) was exploited in Murakawa [11] to design a convergent linear finite-volume scheme, which avoids fully implicit approximations. A semi-implicit finite-volume discretization for a biofilm model with a nonlocal time integrator was proposed in Rahman and Eberl [12]. Finite-volume schemes for cross-diffusion systems with nonlocal (in space) terms were also analyzed; see, for instance, Anaya [13] for a food chain model and [14] for an epidemic model. Moreover, a finite-volume scheme for a Keller-Segel system with additional cross diffusion and discrete entropy dissipation property was investigated in Bessemoulin-Chatard and Jüngel [15]. All these models, however, do not include volume filling and do not possess the degenerate structure explained before.
The paper is organized as follows. The numerical scheme and the main results are presented in Section 2. In Section 3, the existence and uniqueness of bounded discrete solutions are shown. We prove the discrete entropy inequality and further a priori estimates in Section 4, while Section 5 is concerned with the convergence of the numerical scheme. Numerical experiments are given in Section 6 in order to illustrate the order of convergence and the long time behavior of the scheme. For the compactness arguments, we need two discrete Aubin-Lions lemmas which are proved in the appendix.

2
NUMERICAL SCHEME AND MAIN RESULTS

Notations and definitions
We summarize our general hypotheses on the data: H3 Background charge: f ∈ L ∞ (Ω). H4 Initial and boundary data: Remark 1 (Discussion of the assumptions). Assumption (A1) is supposed for simplicity only. Mixed Dirichlet-Neumann boundary conditions can be included in the analysis (see, e.g., [6]), but the proofs become even more technical. Mixed boundary conditions are chosen in the numerical experiments; therefore, the numerical scheme is defined for that case. Assumption (A2) is needed for the derivation of an upper bound for the solvent concentration. Indeed, when D i = D for all i, summing (1) over i = 1, …, n gives On the discrete level, we replace u 0 w Φ by an upwind approximation. This allows us to apply the discrete maximum principle showing that u 0 ≥ 0 and hence u = (u 1 , … , u n ) ∈  with  defined in (6). Finally, Assumption (A3) is needed to derive a discrete version of the entropy inequality. Without the drift terms, the upwinding value does not depend on the index of the species, which simplifies some expressions; see Remark 2.
For the definition of the numerical scheme for (1) and (2), we need to introduce a suitable discretization of the domain Ω and the interval (0, T). For simplicity, we consider a uniform time discretization with time step △ t > 0, and we set t k = k△ t for k = 1, …, N, where T > 0, N ∈ N are given and △ t = T/N. The domain Ω is discretized by a regular and admissible triangulation in the sense of [16,Defin. 9.1]. The triangulation consists of a family  of open polygonal convex subsets of Ω (so-called cells), a family  of edges (or faces in three dimensions), and a family of points (x K ) K∈ associated to the cells. The admissibility assumption implies that the straight line between two centers of neighboring cells x K x L is orthogonal to the edge = K|L between two cells K and L. The condition is satisfied by, for instance, triangular meshes whose triangles have angles smaller than /2 [[16], Exam. 9.1] or Voronoi meshes [16,Exam. 9.2].
We assume that the family of edges  can be split into internal and external edges  =  int ∪  int with  int = { ∈  : ⊂ Ω} and  ext = { ∈  : ⊂ Ω}. Each exterior edge is assumed to be an element of either the Dirichlet or Neumann boundary, that is.  ext =  D ext ∪  N ext . For given K ∈  , we define the set  K of the edges of K, which is the union of internal edges and edges on the Dirichlet or Neumann boundary, and we set The size of the mesh is defined by h( ) = sup{diam(K) ∶ K ∈  }. For ∈  int with = K|L, we denote by d = d(x K , x L ) the Euclidean distance between x K and x L , while for ∈  ext , we set d = d(x K , ). For a given edge ∈ , the transmissibility coefficient is defined by where m( ) denotes the Lebesgue measure of . We impose a regularity assumption on the mesh: There exists > 0 such that for all K ∈  and ∈  K , it holds that This hypothesis is needed to apply discrete functional inequalities (see [16,17]) and a discrete compactness theorem (see [18]).
It remains to introduce suitable function spaces for the numerical discretization. The space   of piecewise constant functions is defined by The (squared) discrete H 1 norm on this space is given by The discrete H −1 norm is the dual norm with respect to the L 2 scalar product, Then Finally, we introduce the space   ,△ t of piecewise constant in time functions with values in   , equipped with the discrete L 2 (0, T;H 1 (Ω)) norm .
For the numerical scheme, we introduce some further definitions. Let u i ∈   with values u i, on the Dirichlet boundary (i = 1, …, n). Then we introduce The numerical fluxes  K, should be consistent approximations to the exact fluxes through the edges ∫  ⋅ ds. We impose the conservation of the numerical fluxes  K, +  L, = 0 for edges = K|L, requiring that they vanish on the Neumann boundary edges,  K, = 0 for ∈  N ext,K . Then the discrete integration-by-parts formula becomes for u ∈   ∑ When Ω = Γ N , this formula simplifies to ∑

Numerical scheme
We need to approximate the initial, boundary, and given functions on the elements K ∈  and edges ∈ : The numerical scheme is as follows. Let K ∈  , k ∈{1, …, N}, i = 1, …, n, and u k−1 i,K ≥ 0 be given. Then the values u k i,K are determined by the implicit Euler scheme where the fluxes  k i,K, are given by the upwind scheme where is defined in (9), and  k i,K, is the "drift part" of the flux, for i = 1, …, n. Observe that we employed a double upwinding: one related to the electric potential, definingû k 0, ,i , and another one related to the drift part of the flux,  k i,K, . The potential is computed via . (20) We recall that the numerical boundary conditions are given by u i, and Φ for ∈  D ext . We denote by u i, ,△ t , Φ  ,△ t the functions in   ,△ t associated to the values u k i,K and Φ k K , respectively. Moreover, when dealing with a sequence of meshes ( m ) m and a sequence of time steps Remark 2 (Simplified numerical scheme)When Assumptions (A1)-(A3) hold, the numerical scheme simplifies to where u k 0,K and u k 0, are defined in (17), and the definition of u k i, simplifies to In the definition of u k i, , the upwinding value does not depend on i anymore such that This property is needed to control the sum ∑ n i=1 u k i, from below in the proof of the discrete entropy inequality; see (35). Finally, we are able to reformulate the discrete fluxes such that we obtain a discrete version of (8) (without the drift part): This formulation is needed in the convergence analysis.

Main results
Since our scheme is implicit and nonlinear, the existence of an approximate solution is nontrivial. Therefore, our first result concerns the well-posedness of the numerical scheme.
Under Assumption (A2), it follows that , and we can apply the discrete minimum principle, which then implies an L ∞ bound for u k i . This bound allows us to apply a topological degree argument; see [19,20]. For the uniqueness proof, we additionally need Assumption (A3), since we use the entropy method of Gajewski [7], and it seems that this method cannot be applied to cross-diffusion systems including drift terms [8]. The idea is to prove first the uniqueness of u k 0 , which solves a discrete nonlinear equation, and then to show the uniqueness of u k showing that it is monotone in k, such that a discrete Gronwall argument implies that u k = v k .
The second result shows that the scheme preserves a discrete version of the entropy inequality.
Assumption (A3) is required to estimate the expression ∑ n i=1 u k i, . In the continuous case, this sum equals 1 − u 0 . On the discrete level, this identity cannot be expected since the value of u k i, depends on the upwinding value; see (18). If the drift part vanishes, the upwinding value does not depend on i, as mentioned in remark 2, and we can derive the estimate Note that the entropy production I k is the discrete counterpart of (7).
The main result of this paper is the convergence of the approximate solutions to a solution to the continuous cross-diffusion system.
where u is a weak solution to (1), (3)-(5) (with Γ N = Ω), that is, for all ∈ C ∞ 0 (Ω×[0, T)) and i = 1, …, n, The compactness of the concentrations follows from the discrete gradient estimates derived from the entropy inequality (25), for which we need Assumption (A3). By the discrete Aubin-Lions lemma [21], we conclude the strong convergence of the sequence (u 1∕2 0,m ). The difficult part is to show the strong convergence of (u 1∕2 0,m u i,m ), since there is no control on the discrete gradient of u i, m . The idea is to apply a discrete Aubin-Lions lemma of "degenerate" type, proved in lemma 10 in the appendix.

L ∞ bounds and existence of solutions
In order to prove the existence of solutions to (15)-(20), we first consider a truncated problem. This means that we truncate the expressions in (18); more precisely, we consider Scheme (15), (16), and (20) with where z + = max {0, z} for z ∈ R and i = 1, …, n. We show that this truncation is, in fact, not needed if the initial data are nonnegative. In the following let (H1)-(H4) hold.

Proof.
We proceed by induction. For k = 0, the nonnegativity holds because of our assumptions on the initial data. Assume that u k−1 i,L ≥ 0 for all L ∈  . Then let u k i,K = min{u k i,L ∶ L ∈  } for some K ∈  and assume that u k i,K < 0. The scheme writes as By assumption, Hence, the right-hand side of (29) nonnegative. However, the left-hand side is negative, which is a contradiction. We infer that u k i,K ≥ 0 and consequently, u k i,L ≥ 0 for all L ∈  . When the initial data are positive, similar arguments show the positivity of u k i,L for L ∈  . ▪ We are able to show the nonnegativity of u k 0,K = 1 − ∑ n i=1 u k i,K only if the diffusion coefficients are the same. The reason is that we derive an equation for u k 0,K by summing (15) for i = 1, …, n, and this gives an equation for u k 0,K only if D i = D for all i = 1, …, n.
Proof. Again, we proceed by induction. The case k = 0 follows from the assumptions. Assume that u k−1 0,L ≥ 0 for all L ∈  . Then let u k 0,K = min{u k 0,L ∶ L ∈  } for some K ∈  and assume that u k 0,K < 0. Summing Equations (15) from i = 1, …, n, we obtain since u k 0, ≥ 0 and u k i, ≥ 0 by construction and D K, (u k 0 ) ≥ 0 because of the minimality property of u k 0,K . The remaining expression is nonnegative: However, the left-hand side of (30) is negative, by induction hypothesis, which gives a contradiction. ▪ Lemmas 4 and 5 imply that we may remove the truncation in (28). Moreover, by definition, we if the initial and boundary data are positive, u k K ∈ .
Proof. We argue by induction. For k = 0, we have u 0 K ∈  by assumption. The function Φ 0 is uniquely determined by scheme (20), as this is a linear system of equations with positive definite matrix. Assume the existence of a solution (u k − 1 , Φ k − 1 ) with u k−1 K ∈ . Let m ∈ N be the product of the number of species n and the number of cells K ∈  . For given K ∈  and i = 1, …, n, we define the function where u 0, K , u i, , u 0, , andû 0, i are defined in (28), and Φ is uniquely determined by (20).

Uniqueness of solutions
The proof of Theorem 1 is completed when we show the uniqueness of solutions to scheme (15)- (20) under the additional conditions (A1) and (A3). Recall that in this case, the scheme is given by (21)-(22), Step 1: uniqueness for u 0 . If k = 0, the solution is uniquely determined by the initial condition. Assume that u k−1 0 is given. Thanks to Assumptions (A2) and (A3), the sum of (21) and (22) for i = 1, …, n gives an equation for u k 0 = 1 − ∑ n i=1 u k i (in the following, we omit the superindices k): where we used (23) We multiply this equation by w 0, K /D, sum over K ∈  , and use discrete integration by parts (14): The first two terms on the right-hand side are clearly nonnegative. We infer from the elementary inequality (y|y| − z|z|)(y − z) ≥ 0 for y, z ∈ R, which is a consequence of the monotonicity of z  → z|z|, that the third term is nonnegative, too. Consequently, the three terms must vanish and this implies that w 0, K = 0 for all K ∈  . This shows the uniqueness for u 0 .
Step 2: uniqueness for u i . Let u 0 be the uniquely determined solution from the previous step and let u k = (u k 1 , … , u k n ) and v k = (v k 1 , … , v k n ) be two solutions to (15). Similarly as in Gajewski [7], we introduce the semimetric and h (z) = (z + )(log(z + ) − 1) + 1. The parameter > 0 is needed since u k i,K or v k i,K may vanish and then the logarithm of u k i,K or v k i,K may be undefined. The objective is to verify that lim →0 (u k , v k ) = 0 by estimating the discrete time derivative of the semimetric, implying that u k = v k . First, we write The function H 1 is convex since ) .
Therefore, a Taylor expansion of ) .
We insert the scheme (21) and (22) and use discrete integration by parts: We claim that S k 1 ≤ 0 and S k 2 ≤ 0. Indeed, with the definition H 2 (a, b) = (a −b)(log(a + )−log(b + )), we can reformulate S k 1 as The Hessian of H 2 , is positive semidefinite. Therefore, performing a Taylor expansion up to second order, we see that S k 1 ≤ 0. Next, we show that S k 2 ≤ 0. For this, we assume without loss of generality for some fixed = K|L that u k 0,K ≤ u k 0,L . By definition of the scheme, u k i, = u k i,K and v k i, = v k i,K . Set H 3 (a, b) = (a + )(log(a + ) − log(b + )). The term in the curly bracket in S k 2 then takes the form The Hessian of H 3 , is also positive semidefinite, showing that (31) is nonpositive as u k 0,K − u k 0,L ≤ 0. If u k 0,K > u k 0,L , both factors of the product (31) change their sign, so that we arrive at the same conclusion. Hence, S k 2 ≤ 0. We conclude that Since d (u 0 , v 0 ) = 0, we find after resolving the recursion that As the densities u i,K are nonnegative and bounded by 1 for all K ∈  , for all ≥ 0 and for all 1 ≤ i ≤ n, it is clear that ∑ k =1 S 3 → 0 when → 0. Then, we may perform the limit → 0 in the previous inequality yielding d (u

. A Taylor expansion as in Zamponi and Jüngel [[8], end of Section 6] shows that
We infer that u k = v k , finishing the proof.

Proof of Theorem 2
The idea is to multiply (15) by log(u k, i,K ∕u k, 0,K ), where u k, i,K ≔ u k i,K + for i = 0, …, n. The regularization is necessary to avoid issues when the concentrations vanish. After this multiplication, we sum the equations over i = 1, …, n and K ∈  and use discrete integration by parts to obtain where The convexity of h(z) = z(logz − 1) In order to estimate the remaining terms, we recall two elementary inequalities. Let y, z > 0. Then, by the Cauchy-Schwarz inequality, and by the concavity of the logarithm, Inequality (33) shows that We use the definition of u k We rewrite B 1 by using the abbreviation u k, i, = u k i, + : We apply inequality (34) to B 11 . Indeed, if u k 0,K ≤ u k 0,L , we have u k i, = u k i,K and we use the first inequality in (34). If u k 0,K > u k 0,L then u k i, = u k i,L and we employ the second inequality in (34). In both cases, it follows that Finally, we consider B 2 . In view of Assumption (A3), Equation (23) gives and therefore, by (33), The last expression cancels with A 2 such that Putting together the estimates for A 0 , A 1 , B 1 , and A 2 + B 2 , we deduce from (32) that Since the right-hand side converges to zero as → 0, we infer that (25) holds.

A priori estimates
For the proof of the convergence result, we need estimates uniform in the mesh size h( ) and time step △ t. The scheme provides uniform L ∞ bounds. Further bounds are derived from the discrete entropy inequality of theorem 2. We introduce the discrete time derivative for functions v ∈   ,△ t by
Proof. We claim that estimates (37) follow from the discrete entropy inequality (25). Indeed, we sum (25) over k = 1, …, N to obtain Since the entropy at time t = 0 is bounded independently of the discretization, we infer immediately the bound for u 1∕2 0 in   ,△ t . For the bound on u 1∕2 0 u i in   ,△ t , we observe that Therefore, together with the L ∞ bounds on Then, summing over k = 0, …, N and using the estimates from the entropy inequality, we achieve the bound on u 1∕2 0 u i . It remains to prove estimate (38). To this end, let ∈   be such that ‖ ‖ 1, = 1 and let k ∈ {1, …, N} and i ∈ {1, …, n}. We multiply the Scheme (21) by Φ K and we sum over K ∈  . Using successively discrete integration by parts, the rewriting of the numerical fluxes (24), the Cauchy-Schwarz inequality, and the L ∞ bounds on u i , we compute . This shows that, for i = 1, …, n, as a consequence of (37). The estimate for follows from those for i = 1, …, n, completing the proof. ▪

CONVERGENCE OF THE SCHEME
In this section, we establish the convergence of the sequence of approximate solutions, constructed in theorem 1, to a weak solution to (1), that is, we prove theorem 3.

Compactness of the approximate solutions
In order to achieve the convergence in the fluxes, we proceed as in Chainais-Hillairet and coworkers [22] by defining the approximate gradient on a dual mesh. For = K | L ∈  int , we define the new cell T KL as the cell with the vertexes x K , x L and those of . For ∈  ext ∩  K , we define T K as the cell with vertex x K and those of . Then Ω can be decomposed as where  K denotes the set of neighboring cells of K.
where n KL denotes the unit normal on = K|L oriented from K to L. To simplify the notation, we set m ≔  m ,△ t m . The solution to the approximate scheme (21) and (22) is called u 0, m , u 1, m , …, u n, m .
Proof. First, we claim that (u 0, m ) is uniformly bounded in   ,△ t . Indeed, by the L ∞ bounds and estimate (37), By estimate (38), ( △ t t u 0,m ) is uniformly bounded. Therefore, by the discrete Aubin-Lions lemma (see lemma 9 in the appendix), we conclude the existence of a subsequence (not relabeled) such that the first convergence in (40) holds. The strong convergence implies (up to a subsequence) that u 0, m → u 0 pointwise in Ω T and consequently u Applying the discrete Aubin-Lions lemma of "degenerate" type (see lemma 10 in the appendix) to y m = u 1∕2 0,m and z m = u i, m for fixed i ∈ {1, …, n}, we deduce convergence (42). Finally, convergence (43) is a consequence of (42) and the weak compactness of (u 1∕2 0,m u i,m ), thanks to the uniform bound in (37).

5.2
The limit m → ∞ We finish the proof of theorem 3 by verifying that the limit function u = (u 1 , …, u n ), as defined in lemma 8, is a weak solution in the sense of the theorem. Let ∈ C ∞ 0 (Ω × [0, T)) and let m ∈ N be large enough such that supp ⊂ Ω × [0, (N m − 1)△ t m ) (recall that T = N m △ t m ). For the limit, we follow the strategy used, for instance, in Chainais-Hillairet and coworkers [22] and introduce the following notations: The convergence results of lemma 8 show that, as m → ∞, Next, setting k K = (x K , t k ), we multiply Scheme (21) by △ t m k−1 K and sum over K ∈  m and k = 1, …, N m . Then where, omitting the subscript m from now on to simplify the notation, The aim is to show that F i0 (m) − F i (m) → 0 as m → ∞ for i = 1, 2, 3. Then, because of (46), F 10 (m) + DF 20 (m) − DF 30 (m) → 0, which finishes the proof. We start by verifying that F 10 (m) − F 1 (m) → 0. For this, we rewrite F 1 (m) and F 10 (m), using N K = 0: In view of the regularity of and the uniform L ∞ bound on u i , we find that Using discrete integration by parts, the second integral becomes where we have decomposed (u k 0, ) 1∕2 = (u k 0,K ) 1∕2 + ((u k 0, ) 1∕2 − (u k 0,K ) 1∕2 ), i.e.
Furthermore, we write F 20 (m) = G 1 (m) + G 2 (m), where The aim is to show that F 21 (m) − G 1 (m) → 0, F 22 (m) → 0, and G 2 (m) → 0. This implies that First we notice that, due to the admissibility of the mesh and the regularity of , by taking the mean value over T KL , where the constant C > 0 only depends on . It yields where the last estimate follows from the Cauchy-Schwarz inequality. This proves that  It remains to analyze the expressions F 22 (m) and G 2 (m). To this end, we remark that d ≤ h( ) and hence, together with the regularity of , and the Cauchy-Schwarz inequality, The term G 2 (m) can be estimated in a similar way. Finally, we need to show that |F 30 (m) − F 3 (m)| → 0. The is completely analogous to the previous arguments, since Summarizing, we have proved that |F i0 (m) − F i (m)| → 0 for i = 1, 2, 3, and since F 1 (m) + DF 2 (m) − DF 3 (m) = 0, the convergence (45) shows that u solves (27). This completes the proof of Theorem 3.

NUMERICAL EXPERIMENTS
We present numerical simulations of a calcium-selective ion channel in two space dimensions to illustrate the dynamical behavior of the ion transport model. Numerical simulations in one space dimension can be found in Burger and coworkers [23] for stationary solutions and in Gerstenmayer and Jüngel [6] for transient solutions. The channel is modeled as in Gillespie coworkers [24]. The selectivity of the channel is obtained by placing some confined oxygen ions (O 1/2− ) inside the channel region. These ions contribute to the permanent charge density f = −u ox /2 in the Poisson equation, but also to the total sum of the concentrations. We consider three further types of ions: calcium (Ca 2+ , u 1 ), sodium (Na + , u 2 ), and chloride (Cl − , u 3 ). While the concentrations of these ion species satisfy the evolution equations (1), the oxygen concentration is constant in time and given by the piecewise linear function The simulations are performed with the full set of Equations (1)-(2) without assuming (A1)-(A3). The finite-volume Scheme (15)- (20) is implemented using MATLAB, version R2015a. The nonlinear system defined by the implicit scheme is solved with a full Newton method in the variables u 0 , u 1 , u 2 , u 3 , Φ for every time step. The computations are done with a fixed time step size △ t = 10 −3 until the stationary state is approximately reached, that is, until the discrete L 2 norm between the solutions at two consecutive time steps is smaller than 10 −12 . We employ an admissible mesh with 4,736 elements generated by the MATLAB command initmesh, which produces Delauney meshes. As initial data, piecewise linear functions that connect the boundary values are chosen for the ion concentrations, while the initial potential is computed from the Poisson equation using the initial concentrations as charge density. Figures 2 and 3 show the concentration profiles and the electric potential after 50 and 1,400 time steps, respectively. The equilibrium is approximately reached after 1,653 time steps. The profiles depicted in Figure 3 are already very close to the stationary state and correspond qualitatively well to the one-dimensional stationary profiles presented in Burger and coworkers [23]. We observe that during the evolution, sodium inside the channel is replaced by the stronger positively charged calcium ions. For higher initial calcium concentrations, the calcium selectivity of the channel acts immediately. The piecewise linear oxygen concentration used in our simulations is clearly not optimal for describing the complex structure of a selectivity filter. As an attempt for a more realistic model, we represent the eight confined oxygen ions as fixed circular discs with a scaled radius of 0.03 (corresponding to 15 nm) and centers inside the channel region. On these discs, the oxygen concentration is set to u ox, max , otherwise it equals zero. Figure 4 shows the concentration profiles and the potential in the equilibrium state. Again, we observe that the calcium ions are selected over the sodium ions, while chloride is removed from the channel.
The simulations suggest that the solution tends towards a steady state as t → ∞. The large-time behavior can be quantified by computing the relative entropy E k with respect to the stationary solution, where and (u ∞ i,K , Φ ∞ ) is the steady state determined from the boundary data. Figure 5 shows that the relative entropy as well as the discrete L 1 norms of the concentrations and electric potential decay with exponential rate. Interestingly, after some initial phase, the convergence is rather slow and increases after this intermediate phase. This phase can be explained by the degeneracy at u 0 = 0, which causes a small entropy production slowing down diffusion. Indeed, as shown in Gerstenmayer and Jüngel [6] for the one-dimensional setting, a small change in the oxygen concentration may prolong the intermediate phase of slow convergence drastically.
The scaled permittivity constant in the Poisson equation has the value 2 = 4.68 × 10 −4 . Therefore, the drift term is moderately convective (in the sense that the modulus of the electric field | Φ| is moderately large). Figure 6 shows the entropy decay and the L 1 error for a smaller value of = 2 × 10 −4 . It turns out that the scheme is still entropy dissipative, but the entropy decay is much slower and the L 1 error decreases with smaller rate.
Since Assumptions (A1)-(A3) are not satisfied in our test case, the convergence result of theorem 3 cannot be applied here. However, we still observe convergence of the numerical solutions. As the exact solution is not known explicitly, we compute a reference solution on a very fine mesh with 75,776 elements and mesh size h( ) ≈ 0.01. This mesh is obtained from the coarse mesh by a regular refinement, dividing the triangles into four triangles of the same shape. The reference solution is compared to approximate solutions on coarser nested meshes. In Figure 7, the errors in the discrete L 1 norm between the reference solution and the solutions on the coarser meshes at two fixed time steps k = 50 and k = 1,400 are plotted. We clearly observe the expected first-order convergence in space.
is the discrete gradient defined in (39), and △ t t is the discrete time derivative defined in (36).
Lemma 9 (Discrete Aubin-Lions). Let ‖⋅‖ 1, m be the norm on   m defined in (11) with the dual norm ‖⋅‖ −1, m given by (12), and let (u m ) ⊂   m ,△ t m be a sequence of piecewise constants in time functions with values in   m satisfying where C > 0 is independent of the size of the mesh and the time step size. Then there exists a subsequence, which is not relabeled, such that, as m → ∞,

Proof.
The result is a consequence of Theorem 3.4 in [18]. To apply this theorem, we have to show that the discrete norms ‖⋅‖ 1, m and ‖⋅‖ −1, m satisfy the assumptions of Lemma 3.1 in [18]: there exists v ∈ L 2 (Ω) such that, up to a subsequence, v m → v in L 2 (Ω). 2 If v m → v strongly in L 2 (Ω) and ‖v m ‖ −1, m → 0 as m → ∞, then v = 0.
The first property is proved in, for instance, [ [21], Lem. 5.6]. Here, we need assumption (10) on the mesh. The second property can be replaced, according to [ [18], Rem. 6], by the condition that ‖⋅‖ 1, m and ‖⋅‖ −1, m are dual norms with respect to the L 2 (Ω) norm, which is the case here. We infer that there exists a subsequence of (u m ), which is not relabeled, such that u m → u strongly in L 2 (Ω T ). The weak convergence of the discrete gradients can be proved as in Lemma 4.4 in Chainais-Hillairet and coworkers [22]. Indeed, the boundedness of ( m u m ) in L 2 implies the convergence to some function ∈ L 2 (Ω T ) (up to a subsequence). In order to show that = u, it remains to verify that for all test functions ∈ C ∞ 0 (Ω T ; R ), This limit follows from the definition of m u m and the regularity of the mesh. We refer to [ [22], Lem. 4.4] for details. ▪ Lemma 10 (Discrete Aubin-Lions of "degenerate" type). Let (y m ) and (z m ) be sequences in   m ,△ t m which are bounded in L ∞ (Ω T ) and let (y m ) be relatively compact in L 2 (Ω T ), i.e., up to a subsequence, y m → y strongly in L 2 (Ω T ) and z m ⇀ * z weakly* in L ∞ (Ω T ). Furthermore, suppose that, for some constant C > 0 independent of m, Then there exists a subsequence which is not relabeled such that y m z m → yz strongly in L 2 (Ω T ) as m → ∞. The L ∞ bounds on z m give By assumption, the sequence (y m ) is relatively compact in L 2 (Ω T ). Therefore, we can apply the inverse of the Kolmogorov-Riesz theorem [ [25], Exerc. 4.34] to conclude that I 21 converges to zero as → 0 uniformly in m.
The analysis of I 22 is more involved. We split the integral in several parts: where the last inequality follows from [ [28], Lems. 4.1 and 4.2]. Let us remark that, for all = K | L ∈  int , we can rewrite (y K ) 2 z K − (y L ) 2 z L = y K + y L 2 (y K z K − y L z L ) + y K z K + y L z L 2 (y K − y L ).
Hence, J 1 ≤ C for some C > 0. An analogous estimation leads to J 2 ≤ C . It remains to estimate the integral J 3 . For this, we use, similar to the treatment of I 21 , the L ∞ bounds on y m and z m : |y m (xt + ) − y m (x, t)|dxdt.
This expression converges to zero uniformly in m because of the relative compactness of (y m ) in L 2 (Ω T ).
We deduce from the previous computations that (48) holds true. Therefore, the product (y m z m ) converges strongly in L 2 (Ω T ), up to some subsequence, and in view of the convergences y m → y strongly in L 2 (Ω T ) and z m ⇀ * z weakly* in L ∞ (0, T; L ∞ (Ω)), the limit of (y m z m ) equals yz, which finishes the proof. ▪