Convergence of a finitevolume scheme for a degenerate crossdiffusion model for ion transport
Abstract.
An implicit Euler finitevolume scheme for a degenerate crossdiffusion system describing the ion transport through biological membranes is analyzed. 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 crossdiffusion system possesses a formal gradientflow structure revealing nonstandard degeneracies, which lead to considerable mathematical difficulties. The finitevolume scheme is based on twopoint flux approximations with “double” upwind mobilities. It preserves the structure of the continuous model like nonnegativity, upper bounds, and entropy dissipation. The degeneracy is overcome by proving a new discrete AubinLions lemma of “degenerate” type. Under suitable assumptions, the existence and uniqueness of bounded discrete solutions, a discrete entropy inequality, and the convergence of the scheme is proved. Numerical simulations of a calciumselective ion channel in two space dimensions indicate that the numerical scheme is of first order.
Key words and phrases:
Ion transport, finitevolume method, gradient flow, entropy method, existence of discrete solutions, convergence of the scheme, calciumselective ion channel.2000 Mathematics Subject Classification:
65M08, 65M12, 35K51, 35K65, 35Q92.1. 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 PoissonNernstPlanck equations [25], 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 onlattice model in the diffusion limit another set of differential equations with fluxes depending on the gradients of all species [10, 27]. These nonlinear crossdiffusion terms are common in multicomponent systems [23, Chapter 4]. In this paper, we propose an implicit Euler finitevolume discretization of the resulting crossdiffusion system. The scheme is designed in such a way that the nonnegativity and upper bound of the concentrations as well as the entropy dissipation is preserved on the discrete level.
More specifically, the evolution of the concentrations and fluxes of the th ion species is governed by the equations
(1) 
for , where is the concentration (volume fraction) of the electroneutral solvent, is a diffusion coefficient, is the (scaled) inverse thermal voltage, and the charge of the th 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 . The electric potential is determined by the Poisson equation
(2) 
where is the (scaled) permittivity constant and is a permanent background charge density. Equations (1) and (2) are solved in a bounded domain .
In order to match experimental conditions, the boundary is supposed to consist of two parts, the insulating part , on which noflux boundary conditions are prescribed, and the union of boundary contacts with external reservoirs, on which the concentrations are fixed. The electric potential is prescribed at the electrodes on . This leads to the mixed DirichletNeumann boundary conditions
(3)  
(4) 
where the boundary data and can be defined on the whole domain . Finally, we prescribe the initial conditions
(5) 
The main mathematical difficulties of equations (1) are the strong coupling and the fact that the diffusion matrix , defined by for and is not symmetric and not positive definite. It was shown in [10, 22] that system (1) possesses a formal gradientflow structure. This means that there exists a (relative) entropy functional with the entropy density
where and , such that (1) can be formally written as
where , for provide a diagonal positive definite matrix, and are the entropy variables, defined by
We refer to [20, Lemma 7] for the computation of .
The entropy structure of (1) is useful for two reasons. First, it leads to bounds for the concentrations. Indeed, the transformation to entropy variables can be inverted, giving with
Then is positive and bounded from above, i.e.
(6) 
This yields bounds without the use of a maximum principle. Second, the entropy structure leads to gradient estimates via the entropy inequality
where the constant depends on the Dirichlet boundary data. Because of
(7) 
we achieve gradient estimates for and . Since may vanish locally, this does not give gradient bounds for , which expresses the degenerate nature of the crossdiffusion system. As a consequence, the flux has to be formulated in the terms of gradients of and only, namely
(8) 
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. This suggests to use the entropy variables as the unknowns, as it was done in our previous work [20] 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 for the numerical discretization. Interestingly, we are still able to prove that the scheme preserves the nonnegativity, upper bound, and entropy inequality. However, the upper bound comes at a price: We need to assume that all diffusion coefficients are the same. Under this assumption, solves a driftdiffusion equation for which the (discrete) maximum principle can be applied. It is not surprising that the bound can be shown only under an additional condition, since crossdiffusion systems usually do not allow for a maximum principle.
The key observation for the numerical discretization is that the fluxes can be written on each cell in a “double” driftdiffusion form, i.e., both and have the structure , where is the diffusion term and is the drift term. We discretize and by using a twopoint flux aproximation with “double” upwind mobilities.
Our analytical results are stated and proved for noflux boundary conditions on . Mixed DirichletNeumann boundary conditions could be prescribed as well, but the proofs would become even more technical. The main results are as follows.

We prove the existence of solutions to the fully discrete numerical scheme (Theorem 1). If 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 [18], recently extended to crossdiffusion systems [20, 28].

Thanks to the “double” upwind structure, the scheme preserves the nonnegativity and upper bound for the concentrations (at least if for all ). Moreover, convexity arguments show that the discrete entropy is dissipated with a discrete entropy production analogous to (7) (Theorem 2). The proof of the discrete entropy inequality only works if the drift term vanishes, since we need to control a discrete version of the sum from below; see the discussion after Theorem 2.

The discrete solutions converge to the continuous solutions 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 AubinLions lemma, which takes into account the nonstandard degeneracy of the equations; see Lemma 10 in the appendix.

Numerical experiments for a calciumselective ion channel in two space dimensions show the dynamical behavior of the solutions and their largetime asymptotics to the equilibrium. The tests indicate that the order of convergence in the norm is one.
In the literature, there exist some results on finitevolume schemes for crossdiffusion systems. An upwind twopoint flux approximation similar to our discretization was recently used in [1] for a seawater intrusion crossdiffusion model. A twopoint flux approximation with a nonlinear positivitypreserving approximation of the crossdiffusion coefficients, modeling the segregation of a twospecies population, was suggested in [4], assuming positive definiteness of the diffusion matrix. The Laplacian structure of the population model (still for positive definite matrices) was exploited in [24] to design a convergent linear finitevolume scheme, which avoids fully implicit approximations. A semiimplicit finitevolume discretization for a biofilm model with a nonlocal time integrator was proposed in [26]. Finitevolume schemes for crossdiffusion systems with nonlocal (in space) terms were also analyzed; see, for instance, [3] for a food chain model and [2] for an epidemic model. Moreover, a finitevolume scheme for a KellerSegel system with additional cross diffusion and discrete entropy dissipation property was investigated in [7]. 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 AubinLions lemmas which are proved in the appendix.
2. Numerical scheme and main results
2.1. Notations and definitions
We summarize our general hypotheses on the data:
(A44)
Domain: ( or ) is an open, bounded, polygonal domain with , .
Parameters: , , , and , .
Background charge: .
Initial and boundary data: , satisfy , and , in for , and .
For our main results, we need additional technical assumptions:
(A44)
, i.e., we impose noflux boundary conditions on the whole boundary.
The diffusion constants are equal, for .
The drift terms are set to zero, .
Remark 1 (Discussion of the assumptions).
Assumption (A1) is supposed for simplicity only. Mixed DirichletNeumann boundary conditions can be included in the analysis (see, e.g., [20]), 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 for all , summing (1) over gives
On the discrete level, we replace by an upwind approximation. This allows us to apply the discrete maximum principle showing that and hence 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)(2), we need to introduce a suitable discretization of the domain and the interval . For simplicity, we consider a uniform time discretization with time step , and we set for , where , are given and . The domain is discretized by a regular and admissible triangulation in the sense of [16, Definition 9.1]. The triangulation consists of a family of open polygonal convex subsets of (socalled cells), a family of edges (or faces in three dimensions), and a family of points associated to the cells. The admissibility assumption implies that the straight line between two centers of neighboring cells is orthogonal to the edge between two cells and . The condition is satisfied by, for instance, triangular meshes whose triangles have angles smaller than [16, Examples 9.1] or Voronoi meshes [16, Example 9.2].
We assume that the family of edges can be split into internal and external edges with and . Each exterior edge is assumed to be an element of either the Dirichlet or Neumann boundary, i.e. . For given , we define the set of the edges of , 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 . For with , we denote by the Euclidean distance between and , while for , we set . For a given edge , the transmissibility coefficient is defined by
(9) 
where denotes the Lebesgue measure of .
We impose a regularity assumption on the mesh: There exists such that for all and , it holds that
(10) 
This hypothesis is needed to apply discrete functional inequalities (see [6, 16]) and a discrete compactness theorem (see [19]).
It remains to introduce suitable function spaces for the numerical discretization. The space of piecewise constant functions is defined by
The (squared) discrete norm on this space is given by
(11) 
The discrete norm is the dual norm with respect to the scalar product,
(12) 
Then
Finally, we introduce the space of piecewise constant in time functions with values in ,
equipped with the discrete norm
For the numerical scheme, we introduce some further definitions. Let with values on the Dirichlet boundary (). Then we introduce
(13)  
The numerical fluxes should be consistent approximations to the exact fluxes through the edges . We impose the conservation of the numerical fluxes for edges , requiring that they vanish on the Neumann boundary edges, for . Then the discrete integrationbyparts formula becomes for
When , this formula simplifies to
(14) 
2.2. Numerical scheme
We need to approximate the initial, boundary, and given functions on the elements and edges :
and we set and .
The numerical scheme is as follows. Let , , , and be given. Then the values are determined by the implicit Euler scheme
(15) 
where the fluxes are given by the upwind scheme
(16) 
where is defined in (9),
(17)  
(18) 
and is the “drift part” of the flux,
(19) 
for . Observe that we employed a double upwinding: one related to the electric potential, defining , and another one related to the drift part of the flux, . The potential is computed via
(20) 
We recall that the numerical boundary conditions are given by and for .
We denote by , the functions in associated to the values and , respectively. Moreover, when dealing with a sequence of meshes and a sequence of time steps , we set , .
Remark 2 (Simplified numerical scheme).
When Assumptions (A1)(A3) hold, the numerical scheme simplifies to
(21)  
(22) 
where and are defined in (17), and the definition of simplifies to
In the definition of , the upwinding value does not depend on anymore such that
(23) 
This property is needed to control the sum 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):
(24) 
This formulation is needed in the convergence analysis. ∎
2.3. Main results
Since our scheme is implicit and nonlinear, the existence of an approximate solution is nontrivial. Therefore, our first result concerns the wellposedness of the numerical scheme.
Theorem 1 (Existence and uniqueness of solutions).
Assumption (A2) is needed to show that is nonnegative. Indeed, summing (15) and (16) over , we obtain
Under Assumption (A2), it follows that , and we can apply the discrete minimum principle, which then implies an bound for . This bound allows us to apply a topological degree argument; see [13, 14]. For the uniqueness proof, we additionally need Assumption (A3), since we use the entropy method of Gajewski [18], and it seems that this method cannot be applied to crossdiffusion systems including drift terms [28]. The idea is to prove first the uniqueness of , which solves a discrete nonlinear equation, and then to show the uniqueness of for by introducing a semimetric for two solutions and and showing that it is monotone in , such that a discrete Gronwall argument implies that .
The second result shows that the scheme preserves a discrete version of the entropy inequality.
Theorem 2 (Discrete entropy inequality).
Assumption (A3) is required to estimate the expression . In the continuous case, this sum equals . On the discrete level, this identity cannot be expected since the value of depends on the upwinding value; see (18). If the drift part vanishes, the upwinding value does not depend on , as mentioned in Remark 2, and we can derive the estimate ; see Section 4.1. Note that the entropy production 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 crossdiffusion system.
Theorem 3 (Convergence of the approximate solution).
Let (H1)(H4) and (A1)(A3) hold and let and be sequences of admissible meshes and time steps, respectively, such that and as . Let be the solution to (21)(22) constructed in Theorem 1. Then there exist functions , satisfying ,
where is a weak solution to (1), (3)(5) (with ), i.e., for all and ,
(27) 
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 AubinLions lemma [17], we conclude the strong convergence of the sequence . The difficult part is to show the strong convergence of , since there is no control on the discrete gradient of . The idea is to apply a discrete AubinLions lemma of “degenerate” type, proved in Lemma 10 in the appendix.
3. Existence and uniqueness of approximate solutions
3.1. 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
(28)  
where for and . We show that this truncation is, in fact, not needed if the initial data are nonnegative. In the following let (H1)(H4) hold.
Lemma 4 (Nonnegativity of ).
Proof.
We proceed by induction. For , the nonnegativity holds because of our assumptions on the initial data. Assume that for all . Then let for some and assume that . The scheme writes as
(29) 
By assumption, . If , we have and if , it follows that . Hence, the righthand side of (29) nonnegative. However, the lefthand side is negative, which is a contradiction. We infer that and consequently, for all . When the initial data are positive, similar arguments show the positivity of for . ∎
We are able to show the nonnegativity of only if the diffusion coefficients are the same. The reason is that we derive an equation for by summing (15) for , and this gives an equation for only if for all .
Lemma 5 (Nonnegativity of ).
Proof.
Again, we proceed by induction. The case follows from the assumptions. Assume that for all . Then let for some and assume that . Summing equations (15) from , we obtain
(30) 
since and by construction and because of the minimality property of . The remaining expression is nonnegative:
However, the lefthand side of (30) is negative, by induction hypothesis, which gives a contradiction. ∎