Comparison of a finite-element and finite-volume scheme for a degenerate cross-diffusion system for ion transport

TL;DR

This paper compares finite-element and finite-volume numerical schemes for simulating ion transport in a degenerate cross-diffusion system, emphasizing structure preservation, convergence, and practical performance.

Contribution

It introduces a structure-preserving finite-element scheme and compares it with a finite-volume scheme for ion transport modeling.

Findings

The finite-element scheme preserves nonnegativity, bounds, and entropy dissipation.

Existence and convergence of discrete solutions are established.

Numerical simulations highlight advantages and drawbacks of both schemes.

Abstract

A structure-preserving implicit Euler finite-element scheme for a degenerate cross-diffusion system for ion transport is analyzed. The scheme preserves the nonnegativity and upper bounds of the ion concentrations, the total relative mass, and it dissipates the entropy (or free energy). The existence of discrete solutions to the scheme and their convergence towards a solution to the continuous system is proved. Numerical simulations of two-dimensional ion channels using the finite-element scheme with linear elements and an alternative finite-volume scheme are presented. The advantages and drawbacks of both schemes are discussed in detail.

Tables

Table 1. Table 1. Difference between the finite-volume and finite-element reference solutions after 600 and after 6000 time steps.

	$u_1$	$u_2$	$u_3$	$u_0$	$\Phi$
$L^{\infty }$ norm, $k=600$	2.2405e-02	2.0052e-02	1.0319e-04	1.6695e-02	1.0600e-01
$L^1$ norm, $k=600$	2.2642e-04	3.0275e-04	1.3776e-05	2.5983e-04	5.1029e-03
$L^{\infty }$ norm, $k=6000$	1.0036e-02	2.3619e-03	1.3677e-04	9.1095e-03	9.5080e-02
$L^1$ norm, $k=6000$	1.4161e-04	7.0981e-05	1.5498e-05	1.5615e-04	4.6543e-03

Table 2. Table 2. Average time needed to compute one time step (in seconds). FE = finite-element scheme, FV = finite-volume scheme.

	${𝒯}_1$	${𝒯}_2$	${𝒯}_3$	${𝒯}_4$	${𝒯}_5$
FE	2.4065e-01	7.9982e-01	2.1125e+00	4.9844e+00	17.7788e+00
FV	6.7707e-03	2.2042e-02	3.0532e-01	1.7660e+00	2.2418e+00

Table 3. Table 3. Dimensionless parameters used for the simulation of the bipolar ion channel and values used for the scaling.

Meaning	Value	Unit
Diffusion coefficients $D_1$, $D_2$	1
Effective permittivity ${\lambda }^2$	1.1713
Effective mobility $\beta$	3.8922
Bath concentrations $\overline{u_1}$, $\overline{u_2}$	0.0016
Confined ion concentration $c_{\max }$	0.2971
Typical length $L_{\text{typ}}$	1e-10	m
Typical concentration $u_{\text{typ}}$	3.7037e+28	Nm-3
Typical voltage ${\Phi }_{\text{typ}}$	0.1	V
Typical diffusion $D_{\text{typ}}$	1.3340e-9	m2s-1

Equations

$$\partial_{t}u_{i}+\operatorname{div}{\mathcal{F}}_{i}=0,\quad{\mathcal{F}}_{i}=-D_{i}\big{(}u_{0}\nabla u_{i}-u_{i}\nabla u_{0}+u_{0}u_{i}\beta z_{i}\nabla\Phi\big{)}\quad\mbox{in }\Omega,\ t>0,$$

$$-{\lambda }^2\Delta \Phi =\sum _{i=1}^n{z}_i{u}_i+f\text{mbox}in\Omega ,$$

$${\mathcal{F}}_i\cdot \nu =0\text{mbox}on{\Gamma }_N,$$

$$\nabla \Phi \cdot \nu =0\text{mbox}on{\Gamma }_N,$$

$$H(u)={\int }_{\Omega }h(u)dx,\text{mbox}whereh(u)=\sum _{i=0}^n{\int }_{{\overline{u}}_i}^{u_i}\log \frac{s}{{\overline{u}}_i}ds+\frac{\beta {\lambda }^2}{2}|\nabla (\Phi -\overline{\Phi }){|}^2.$$

$$\frac{\partial h}{\partial u_i}=w_i-{\overline{w}}_i,\text{mbox}where$$

$$w_i=\log \frac{u_i}{u_0}+\beta z_i\Phi ,{\overline{w}}_i=\log \frac{{\overline{u}}_i}{{\overline{u}}_0}+\beta z_i\overline{\Phi },i=1,\dots ,n,$$

$$\partial_{t}u_{i}(w,\Phi)=\operatorname{div}\bigg{(}\sum_{j=1}^{n}B_{ij}(w,\Phi)\nabla w_{j}\bigg{)},$$

$$B_{ij}(w,\Phi )=D_i{u}_0(w,\Phi )u_i(w,\Phi ){\delta }_{ij},i,j=1,\dots ,n,$$

$$u_i(w,\Phi )=\frac{\exp (w_i-\beta z_i\Phi )}{1+{\sum }_{j=1}^n\exp (w_j-\beta z_j\Phi )},i=1,\dots ,n,$$

$$u(w,\Phi)\in\mathcal{D}:=\bigg{\{}u\in(0,1)^{n}:\sum_{i=1}^{n}u_{i}<1\bigg{\}}.$$

$$\frac{dH}{dt}+\frac{1}{2}{\int }_{\Omega }\sum _{i=1}^n{D}_i{u}_0{u}_i|\nabla w_i{|}^{2}dx\le C,$$

$${\rho }_K\le h_K\le \gamma {\rho }_K\text{mbox}forallK\in {\mathcal{T}}_h,$$

$$\mathcal{S}({\mathcal{T}}_h)$$

$${\mathcal{S}}_D({\mathcal{T}}_h)$$

$$\underset{h\to 0}{\lim }\Vert \varphi -I_h\varphi {\Vert }_{H^1(\Omega )}=0\text{for all }\varphi \in H^2(\Omega ).$$

$$\displaystyle\frac{1}{\tau}\int_{\Omega}\big{(}u(w^{k},\Phi^{k})-u(w^{k-1},\Phi^{k-1})\big{)}\cdot\phi\,dx$$

$$+{\int }_{\Omega }\nabla \varphi :B(w^k,{\Phi }^k)\nabla w^{k}dx+\epsilon {\int }_{\Omega }(w^k-{\overline{w}}_h)\cdot \varphi dx=0,$$

$$\displaystyle\lambda^{2}\int_{\Omega}\nabla\Phi^{k}\cdot\nabla\theta dx=\int_{\Omega}\bigg{(}\sum_{i=1}^{n}z_{i}u_{i}(w^{k},\Phi^{k})+f\bigg{)}\theta dx$$

$$\nabla \varphi :B(w^k,{\Phi }^k)\nabla w^k=\sum _{i=1}^n{D}_i{u}_i(w^k,{\Phi }^k)u_0(w^k,{\Phi }^k)\nabla {\varphi }_i\cdot \nabla w_i^k.$$

$$H(u^k)+\tau {\int }_{\Omega }\nabla (w^k-{\overline{w}}_h):B(w^k,{\Phi }^k)\nabla w^{k}dx+\epsilon \tau \Vert w^k-{\overline{w}}_h{\Vert }_{L^2(\Omega )}^2\le H(u^{k-1}),$$

$$\Vert {\Phi }^k{\Vert }_{H^1(\Omega )}\le C(1+\Vert {\overline{\Phi }}_h{\Vert }_{H^1(\Omega )})$$

$$a(v,\varphi )=F(\varphi )\text{mbox}forall\varphi \in {\mathcal{S}}_D({\mathcal{T}}_h{)}^n,$$

$$a(v,\varphi )$$

$$F(\varphi )$$

$$-\delta {\int }_{\Omega }\nabla \varphi :B(y+{\overline{w}}_h,{\Phi }^k)\nabla {\overline{w}}_{h}dx.$$

$$a(v,v)\ge \epsilon \Vert v{\Vert }_{L^2(\Omega )}^2\ge \epsilon C\Vert v{\Vert }_{H^1(\Omega )}^2.$$

$$\epsilon C\Vert v{\Vert }_{H^1(\Omega )}^2\le a(v,v)=F(v)\le C(\tau )\Vert v{\Vert }_{H^1(\Omega )}$$

$${\int }_{\Omega }(u^k-u^{k-1})\cdot (w^k-{\overline{w}}_h)dx={\int }_{\Omega }(u^k-u^{k-1})\cdot \nabla h(u^k)dx\ge H(u^k)-H(u^{k-1}),$$

$$\Vert u_i^k{\Vert }_{L^{\infty }(\Omega )}+\epsilon \tau \sum _{j=1}^k\Vert w_i^j-{\overline{w}}_{i,h}{\Vert }_{L^2(\Omega )}^2$$

$$\displaystyle\tau\sum_{j=1}^{k}\Big{(}\|(u_{0}^{j})^{1/2}\|_{H^{1}(\Omega)}^{2}+\|u_{0}^{j}\|_{H^{1}(\Omega)}^{2}+\|(u_{0}^{j})^{1/2}\nabla(u_{i}^{j})^{1/2}\|_{L^{2}(\Omega)}^{2}\Big{)}$$

$$\nabla (w^k-{\overline{w}}_h):B(w^k,{\Phi }^k)\nabla w^k$$

$$-\frac{D_{\max }}{2}\sum _{i=1}^n|\nabla {\overline{w}}_i{|}^2,$$

$$H(u^k)$$

$$\le H(u^{k-1})+\tau \frac{D_{\min }}{2}\sum _{i=1}^n|\beta z_i\nabla {\Phi }^k{|}^{2}dx+\tau \frac{D_{\max }}{2}{\int }_{\Omega }\sum _{i=1}^n|\nabla {\overline{w}}_{i,h}{|}^{2}dx.$$

Full text

The authors acknowledge partial support from the Austrian Science Fund (FWF), grants F65, P27352, P30000, and W1245

1. Introduction

Ion channels are pore-forming proteins that create a pathway for charged ions to pass through the cell membrane. They are of great biological importance since they contribute to processes in the nervous system, the coordination of muscle contraction, and the regulation of secretion of hormones, for instance. Ion-channel models range from simple systems of differential equations [18] as well as Brownian and Langevin dynamics [20, 26] to the widely used Poisson–Nernst–Planck model [11]. The latter model fails in narrow channels since it neglects the finite size of the ions. Finite-size interactions can be approximately captured by adding suitable chemical potential terms [16, 27], for instance. In this paper, we follow another approach. Starting from a random walk on a lattice, one can derive in the diffusion limit an extended Poisson–Nernst–Planck model, taking into account that ion concentrations might saturate in the narrow channel. This leads to the appearance of cross-diffusion terms in the evolution equations for the ion concentrations [3, 32]. These nonlinear cross-diffusion terms are common in diffusive multicomponent systems [22, Chapter 4]. A lattice-free approach, starting from stochastic Langevin equations, can be found in [2]. The scope of this paper is to present a new finite-element discretization of the degenerate cross-diffusion system and to compare this scheme to a previously proposed finite-volume method [5].

The dynamics of the ion concentrations $u=(u_{1},\ldots,u_{n})$ is governed by the evolution equations

[TABLE]

where $u_{0}=1-\sum_{i=1}^{n}u_{i}$ denotes the solvent concentration, $D_{i}>0$ is the diffusion constant, $z_{i}$ the ion charge, and $\beta$ a mobility parameter. To be precise, $u_{i}$ is the mass fraction of the $i$th ion, and we refer to $\sum_{i=0}^{n}u_{i}=1$ as the total relative mass, just meaning that the ion-solvent mixture is saturated. The electric potential $\Phi$ is self-consistently given by the Poisson equation

[TABLE]

with the permanent charge density $f=f(x)$ and the scaled permittivity constant $\lambda^{2}$. The equations are solved in a bounded domain $\Omega\subset{\mathbb{R}}^{d}$ with smooth boundary $\partial\Omega$. Equations (1) are equipped with initial data $u(0)=u^{I}$ satisfying $0<\sum_{i=1}^{n}u_{i}^{I}<1$. The boundary $\partial\Omega$ consists of an insulating part $\Gamma_{N}$ and the union $\Gamma_{D}$ of contacts with external reservoirs:

[TABLE]

System (1)-(2) can be interpreted as a generalized Poisson–Nernst–Planck model. The usual Poisson–Nernst–Planck equations [11] follow from (1) by setting $u_{0}=\mbox{const}$. In the literature, there are several generalized versions of the standard model. For instance, adding a term involving the relative velocity differences in the entropy production leads to cross-diffusion expressions different from (1) [19]. This model, however, does not take into account effects from the finite ion size. Thermodynamically consistent Nernst–Planck models with cross-diffusion terms were suggested in [9], but the coefficients differ from (1). The model at hand was derived in [3, 32] from a lattice model taking into account finite-size effects.

Model (1)-(4) contains some mathematical difficulties. First, its diffusion matrix $A(u)=(A_{ij}(u))\in{\mathbb{R}}^{n\times n}$, given by $A_{ij}(u)=D_{i}u_{i}$ for $i\neq j$ and $A_{ii}(u)=D_{i}(u_{0}+u_{i})$ for $i=1,\ldots,n$ is generally neither symmetric nor positive definite. Second, it degenerates in regions where the concentrations vanish. Third, the standard maximum principle cannot be applied to achieve $0\leq u_{i}\leq 1$ for $i=1,\ldots,n$. In the following, we explain how these issues can be solved.

The first difficulty can be overcome by introducing so-called entropy variables $w_{i}$ defined from the entropy (or, more precisely, free energy) of the system,

[TABLE]

Indeed, writing equations (1) in terms of the entropy variables $w_{1},\ldots,w_{n}$, given by

[TABLE]

it follows that

[TABLE]

where the new diffusion matrix $B=(B_{ij}(w,\Phi))\in{\mathbb{R}}^{n\times n}$ with

[TABLE]

is symmetric and positive semidefinite (in fact, it is even diagonal). This procedure has a thermodynamical background: The quantities $\partial h/\partial u_{i}$ are known as the chemical potentials, and $B$ is the so-called mobility or Onsager matrix [8].

The transformation to entropy variables also solves the third difficulty. Solving the transformed system (7) for $w=(w_{1},\ldots,w_{n})$, the concentrations are given by

[TABLE]

showing that $u_{i}$ is positive and bounded from above:

[TABLE]

Moreover, the entropy structure leads to gradient estimates via the entropy inequality

[TABLE]

where the constant $C>0$ depends on the Dirichlet boundary data.

Still, we have to deal with the second difficulty, the degeneracy. It is reflected in the entropy inequality since we lose the gradient estimate if $u_{i}=0$ or $u_{0}=0$. This problem is overcome by using the “degenerate” Aubin-Lions lemma of [21, Appendix C] or its discrete version in [5, Lemma 10].

These ideas were employed in [3] for $n=2$ ion species and without electric potential to show the global existence of weak solutions. The existence result was extended to an arbitrary number of species in [21, 33], still excluding the potential. A global existence result for the full problem (1)-(4) was established in [15].

We are interested in devising a numerical scheme which preserves the structure of the continuous system, like nonnegativity, upper bounds, and the entropy structure, on the discrete level. A first result in this direction was presented in [5], analyzing a finite-volume scheme preserving the aforementioned properties. However, the scheme preserves the nonnegativity and upper bounds only if the diffusion coefficients $D_{i}$ are all equal, and the discrete entropy is dissipated only if additionally the potential term vanishes. In this paper, we propose a finite-element scheme for which the structure preservation holds under natural conditions.

Before we proceed, we briefly discuss some related literature. While there are many results for the classical Poisson–Nernst–Planck system, see for example [25, 29], there seems to be no numerical analysis of the ion-transport model (1)-(4) apart from the finite-volume scheme in [5] and simulations of the stationary equations in [4]. Let us mention some other works on finite-element methods for related cross-diffusion models. In [1], a convergent finite-element scheme for a cross-diffusion population model was presented. The approximation is not based on entropy variables, but a regularization of the entropy itself that is used to define a regularized system. The same technique was employed also in [14]. A lumped finite-element method was analyzed in [13] for a reaction-cross-diffusion equation on a stationary surface with positive definite diffusion matrix. In [23], an implicit Euler Galerkin approximation in entropy variables for a Poisson–Maxwell–Stefan system was shown to converge. Recently, an abstract framework for the numerical approximation of evolution problems with entropy structure was presented in [10]. The discretization is based on a discontinuous Galerkin method in time and a Galerkin approximation in space. When applied to cross-diffusion systems, this approach also leads to an approximation in entropy variables that preserves the entropy dissipation. However, neither the existence of discrete solutions nor the convergence of the scheme are discussed.

Our main results are as follows:

• •

We propose an implicit Euler finite-element scheme for (1)-(4) in entropy variables with linear finite elements (Section 2). The scheme preserves the nonnegativity of the concentrations and the upper bounds, the total relative mass, and it dissipates the discrete entropy associated to (5) if the boundary data are in thermal equilibrium; see Fthe Remark 1.

• •

We prove the existence of discrete solutions (Lemma 1) and their convergence to the solution to (1)-(4) when the approximation parameters tend to zero (Theorem 3). The convergence rate can be only computed numerically and is approximately of second order (with respect to the $L^{2}$ norm).

• •

The finite-element scheme and the finite-volume scheme of [5] (recalled in Section 3) are applied to two test cases in two space dimensions: a calcium-selective ion channel and a bipolar ion channel (Section 4). Static current-voltage curves show the rectifying behavior of the bipolar ion channel.

• •

The advantages and drawbacks of both schemes are discussed (Section 5). The finite-element scheme allows for structure-preserving properties under natural assumptions, while the finite-volume scheme can be analyzed only under restrictive conditions. On the other hand, the finite-volume scheme allows for vanishing initial concentrations and faster algorithms compared to the finite-element scheme due to the highly nonlinear structure of the latter formulation.

2. The finite-element scheme

2.1. Notation and assumptions

Before we define the finite-element discretization, we introduce our notation and make precise the conditions assumed throughout this section. We assume:

{labeling}

(A44)

Domain: $\Omega\subset{\mathbb{R}}^{d}$ ($d=2$ or $d=3$) is an open, bounded, polygonal domain with $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}\in C^{0,1}$, $\Gamma_{D}\cap\Gamma_{N}=\emptyset$, $\Gamma_{N}$ is open in $\partial\Omega$, and $\mbox{meas}(\Gamma_{D})>0$.

Parameters: $T>0$, $D_{i}>0$, $\beta>0$, and $z_{i}\in{\mathbb{R}}$, $i=1,\ldots,n$.

Background charge: $f\in L^{\infty}(\Omega)$.

Initial and boundary data: $u_{i}^{I}\in H^{2}(\Omega)$ and $\overline{u}_{i}\in H^{2}(\Omega)$ satisfy $u_{i}^{I}>0$, $\overline{u}_{i}>0$ for $i=1,\ldots,n$, $1-\sum_{i=1}^{n}u_{i}^{I}>0$, $1-\sum_{i=1}^{n}\overline{u}_{i}>0$ in $\Omega$, and $\overline{\Phi}\in H^{2}(\Omega)\cap L^{\infty}(\Omega)$.

The $H^{2}$ regularity of the initial and boundary data ensures that the standard interpolation converges to the given data, see (10) below.

We consider equations (1) on a finite time interval $(0,T)$ with $T>0$. For simplicity, we use a uniform time discretization with time step $\tau>0$ and set $t^{k}=k\tau$ for $k=1,\ldots,N$, where $N\in{\mathbb{N}}$ is given and $\tau=T/N$.

For the space discretization, we introduce a family ${\mathcal{T}}_{h}$ ($h>0$) of triangulations of $\Omega$, consisting of open polygonal convex subsets of $\Omega$ (the so-called cells) such that $\overline{\Omega}=\cup_{K\in{\mathcal{T}}_{h}}\overline{K}$ with maximal diameter $h=\max_{K\in{\mathcal{T}}_{h}}\text{diam}(K)$. We assume that the corresponding family of edges ${\mathcal{E}}$ can be split into internal and external edges, ${\mathcal{E}}={\mathcal{E}}_{\rm int}\cup{\mathcal{E}}_{\rm ext}$ with ${\mathcal{E}}_{\rm int}=\{\sigma\in{\mathcal{E}}:\sigma\subset\Omega\}$ and ${\mathcal{E}}_{\rm ext}=\{\sigma\in{\mathcal{E}}:\sigma\subset\partial\Omega\}$. Each exterior edge is assumed to be an element of either the Dirichlet or Neumann boundary, i.e. ${\mathcal{E}}_{\rm ext}={\mathcal{E}}_{\rm ext}^{D}\cup{\mathcal{E}}_{\rm ext}^{N}$. For given $K\in{\mathcal{T}}_{h}$, we define the set ${\mathcal{E}}_{K}$ of the edges of $K$, which is the union of internal edges and edges on the Dirichlet or Neumann boundary, and we set ${\mathcal{E}}_{K,\rm int}={\mathcal{E}}_{K}\cap{\mathcal{E}}_{\rm int}$.

In the finite-element setting, the triangulation is completed by the set of nodes $\{p_{j}:j\in J\}$. We impose the following regularity assumption on the mesh. There exists a constant $\gamma\geq 1$ such that

[TABLE]

where $\rho_{K}$ is the radius of the incircle and $h_{K}$ is the diameter of $K$.

We associate with ${\mathcal{T}}_{h}$ the usual conforming finite-element spaces

[TABLE]

and $H_{D}^{1}(\Omega)$ is the set of $H^{1}(\Omega)$ functions that vanish on $\Gamma_{D}$ in the weak sense. Let $\{\chi_{j}\}_{j\in J}$ be the standard basis functions for ${\mathcal{S}}({\mathcal{T}}_{h})$ with $\chi_{j}(p_{i})=\delta_{ij}$ for all $i,j\in J$. We define the nodal interpolation operator $I_{h}:C^{0}(\overline{\Omega})\to{\mathcal{S}}({\mathcal{T}}_{h})$ via $(I_{h}v)(p_{j})=v(p_{j})$ for all $v\in{\mathcal{S}}({\mathcal{T}}_{h})$ and $j\in J$. Due to the regularity assumptions on the mesh, $I_{h}$ has the following approximation property (see, e.g., [7, Chapter 3]):

[TABLE]

2.2. Definition of the scheme

To define the finite-element scheme, we need to approximate the initial and boundary data. We set $w^{0}_{i}=I_{h}(\log(u_{i}^{I}/u_{0}^{I}))+\beta z_{i}\Phi^{0}$, where $\Phi^{0}$ is the standard finite-element solution to the linear equation (2) with $u_{i}^{I}$ on the right-hand side. Furthermore, we set $\overline{w}_{h}=I_{h}(\log(\overline{u}_{i}/\overline{u}_{0})+\beta z_{i}\overline{\Phi})$ and $\overline{\Phi}_{h}=I_{h}(\overline{\Phi})$.

The finite-element scheme is now defined as follows. Given $w^{k-1}\in{\mathcal{S}}({\mathcal{T}}_{h})^{n}$ and $\Phi^{k-1}\in{\mathcal{S}}({\mathcal{T}}_{h})$, find $w^{k}-\overline{w}_{h}\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$ and $\Phi^{k}-\overline{\Phi}_{h}\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})$ such that

[TABLE]

for all $\phi\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$ and $\theta\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})$. The symbol “:” signifies the Frobenius matrix product; here, the expression reduces to

[TABLE]

The term involving the parameter $\varepsilon>0$ is only needed to guarantee the coercivity of (11)-(12). Indeed, the diffusion matrix $B(w^{k},\Phi^{k})$ degenerates when $w^{k}_{i}\to-\infty$, and the corresponding bilinear form is only positive semidefinite. To emphasize the dependence on the mesh and $\varepsilon$, we should rather write $w^{(h,\varepsilon,k)}$ instead of $w^{k}$ and similarly for $\Phi^{k}$; however, for the sake of presentation, we will mostly omit the additional superscripts. The original variables are recovered by computing $u^{k}=u(w^{k},\Phi^{k})$ according to (8). Setting $u^{(\tau)}(x,t)=u^{k}(x)$ for $x\in\Omega$, $t\in((k-1)\tau,k\tau]$, $k=1,\ldots,N$, and $u^{(\tau)}(\cdot,0)=I_{h}u^{I}$ as well as similarly for $\Phi^{(\tau)}$, we obtain piecewise constant in time functions.

2.3. Existence of discrete solutions

The first result concerns the existence of solutions to the nonlinear finite-element scheme (11)-(12).

Lemma 1** (Existence of solutions and discrete entropy inequality).**

There exists a solution to scheme (11)-(12) that satisfies the following discrete entropy inequality:

[TABLE]

where $H$ is defined in (5) and $u^{k}=u(w^{k},\Phi^{k})$, $u^{k-1}=u(w^{k-1},\Phi^{k-1})$ are defined in (8).

The proof of the lemma is similar to the proof of Theorem 1 in [15]. The main difference is that in [15], a regularization term of the type $\varepsilon((-\Delta)^{m}w^{k}+w^{k})$ has been added to achieve via $H^{m}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ for $m>d/2$ compactness and $L^{\infty}$ solutions. In the finite-dimensional setting, this embedding is not necessary but we still need the regularization $\varepsilon w^{k}$ to conclude coercivity. We conjecture that this regularization is just technical but currently, we are not able to remove it. Note, however, that we can use arbitrarily small values of $\varepsilon$ in the numerical simulations such that the additional term does not affect the solution practically.

Proof of Lemma 1..

Let $y\in{\mathcal{S}}({\mathcal{T}}_{h})^{n}$ and $\delta\in[0,1]$. There exists a unique solution $\Phi^{k}$ to (12) with $w^{k}$ replaced by $y+\overline{w}_{h}$, satisfying $\Phi^{k}-\overline{\Phi}_{h}\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})$, since the function $\Phi\mapsto u_{i}(y,\Phi)$ is bounded. Moreover, the estimate

[TABLE]

holds for some constant $C>0$.

Next, we consider the linear problem

[TABLE]

where

[TABLE]

The bilinear form $a$ and the linear form $F$ are continuous on ${\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$. The equivalence of all norms on the finite-dimensional space ${\mathcal{S}}_{D}({\mathcal{T}}_{h})$ implies the coercivity of $a$,

[TABLE]

By the Lax-Milgram lemma, there exists a unique solution $v\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$ to this problem. This defines the fixed-point operator $S:{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}\times[0,1]\to{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$, $S(y,\delta)=v$. The inequality

[TABLE]

shows that all elements $v$ are bounded independently of $y$ and $\delta$ and thus, all fixed points $v=S(v,\delta)$ are uniformly bounded. Furthermore, $S(y,0)=0$ for all $y\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$. The continuity of $S$ follows from standard arguments and the compactness comes from the fact that ${\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$ is finite-dimensional. By the Leray-Schauder fixed-point theorem, there exists $v^{k}\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$ such that $S(v^{k},1)=v^{k}$, and $w^{k}:=v^{k}+\overline{w}_{h}$ is a solution to (11).

The discrete entropy inequality (13) is proven by using $\tau(w^{k}-\overline{w}_{h})\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$ as a test function in (11) and exploiting the convexity of $H$,

[TABLE]

which concludes the proof. ∎

Remark 1** (Structure-preservation of the scheme).**

Lemma 1 shows that if the boundary data is in thermal equilibrium, i.e. $\nabla\overline{w}_{h}=0$, then the finite-element scheme (11)-(12) dissipates the entropy (5), i.e. $H(u^{k})\leq H(u^{k-1})$. Moreover, it preserves the invariant region $\mathcal{D}$, i.e. $u^{k}\in\mathcal{D}$, and the mass fraction $u_{i}^{k}$ is nonnegative and bounded by one. The scheme conserves the total relative mass, i.e. $\sum_{i=0}^{n}\|u_{i}^{k}\|_{L^{1}(\Omega)}=1$, which is a direct consequence of the definition of $u_{0}^{k}$. ∎

2.4. Uniform estimates

The next step is the derivation of a priori estimates uniform in the parameters $\varepsilon$, $\tau$, and $h$. To this end, we transform back to the original variable $u^{k}$ and exploit the discrete entropy inequality (13).

Lemma 2** (A priori estimates).**

For the solution to the finite-element scheme from Lemma 1, the following estimates hold:

[TABLE]

for $i=1,\ldots,n$, where here and in the following, $C>0$ is a generic constant independent of $\varepsilon$, $\tau$, and $h$.

Proof.

As the proof is similar to that one in the continuous setting, we give only a sketch. Observe that the definition of the entropy variables implies that $0<u_{i}^{k}<1$ in $\Omega$ for $i=1,\ldots,n$ and $k=1,\ldots,N$. It is shown in the proof of Lemma 6 of [15] that

[TABLE]

where $D_{\rm min}=\min_{i=1,\ldots,n}D_{i}$ and $D_{\rm max}=\max_{i=1,\ldots,n}D_{i}$. Then (13) gives

[TABLE]

We resolve this recursion to find that

[TABLE]

The right-hand side is bounded because of (14), $\tau k\leq T$, and the boundedness of the interpolation operator. Inserting the identity

[TABLE]

the estimates follow. ∎

2.5. Convergence of the scheme

The a priori estimates from the previous lemma allow us to formulate our main result, the convergence of the finite-element solutions to a solution to the continuous model (1)-(4).

Theorem 3** (Convergence of the finite-element solution).**

Let $(u^{(h,\varepsilon,\tau)},\Phi^{(h,\varepsilon,\tau)})$ be an approximate solution constructed from scheme (11)-(12). Set $u^{(h,\varepsilon,\tau)}_{0}=1-\sum_{i}u^{(h,\varepsilon,\tau)}_{i}$. Then there exist functions $u_{0}$, $u=(u_{1},\ldots,u_{n})$, and $\Phi$, satisfying $u(x,t)\in\overline{\mathcal{D}}$ (${\mathcal{D}}$ is defined in (9)), $u_{0}=1-\sum_{i=1}^{n}u_{i}$ in $\Omega$, the regularity

[TABLE]

for $i=1,\ldots,n$, such that as $(h,\varepsilon,\tau)\to 0$,

[TABLE]

and $(u,\Phi)$ are a weak solution to (1)-(4). In particular, for all $\phi\in L^{2}(0,T;H^{1}_{D}(\Omega))$ and $i=1,\ldots,n$, it holds that

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the duality pairing in $H^{1}_{D}(\Omega)^{\prime}$ and $H_{D}^{1}(\Omega)$, and the boundary and initial conditions are satisfied in a weak sense.

Proof.

We pass first to the limit $(\varepsilon,h)\to 0$ and then $\tau\to 0$, since the latter limit can be performed as in the proof of Theorem 1 in [15]. Fix $k\in\{1,\ldots,N\}$ and let $u_{i}^{(\varepsilon,h)}=u_{i}^{(\varepsilon,h,k)}$ and $\Phi^{(\varepsilon,h)}=\Phi^{(\varepsilon,h,k)}$ be the approximate solution from Lemma 1. We set $u_{0}^{(\varepsilon,h)}=1-\sum_{i=1}^{n}u_{i}^{(\varepsilon,h)}$. Using the compact embedding $H^{1}(\Omega)\hookrightarrow L^{2}(\Omega)$ and the a priori estimates from Lemma 2, it follows that there exists a subsequence which is not relabeled such that, as $(\varepsilon,h)\to 0$,

[TABLE]

Combining (17) and (21), we infer that (up to a subsequence)

[TABLE]

Next, let $\phi\in(H^{2}(\Omega)\cap H_{D}^{1}(\Omega))^{n}$. As we cannot use $\phi_{i}$ directly as a test function in (11), we take $I_{h}\phi\in{\mathcal{S}}_{D}({\mathcal{T}}_{h})^{n}$, where $I_{h}$ is the interpolation operator, see (10). In order to pass to the limit in (11), we rewrite the integral involving the diffusion matrix:

[TABLE]

We estimate each of the above summands separately. For the last term, we proceed as follows:

[TABLE]

Similarly as for (24), it follows that

[TABLE]

Then, together with the weak convergence of $\nabla\Phi^{(\varepsilon,h)}$, we infer that

[TABLE]

Since $(u_{i}^{(\varepsilon,h)}u_{0}^{(\varepsilon,h)}\nabla\Phi^{(\varepsilon,h)})$ is bounded in $L^{2}(\Omega)$, this weak convergence also holds in $L^{2}(\Omega)$. Because of the interpolation property (10) and estimate (26),

[TABLE]

Following the arguments of Step 3 in [15, Section 2] (using (24)), we have

[TABLE]

Thus, the limit $(\varepsilon,h)\to 0$ in (25) gives

[TABLE]

Furthermore, we deduce from (23) that

[TABLE]

Thus, passing to the limit $(\varepsilon,h)\to 0$ in scheme (11)-(12) leads to

[TABLE]

for all $\phi_{i}$, $\theta\in H^{2}(\Omega)\cap H_{D}^{1}(\Omega)$. A density argument shows that we can take test functions $\phi_{i}$, $\theta\in H_{D}^{1}(\Omega)$. The a priori estimates from Lemma 2 remain valid in the weak limit.

Now the limit $\tau\to 0$ can be done exactly as in [15, Theorem 1, Step 4], which concludes the proof. ∎

3. The finite-volume scheme

We briefly recall the finite-volume scheme from [5] and summarize the assumptions and results, as this is necessary for the comparison of the finite-element and finite-volume scheme in Section 5.

We assume that Hypotheses (H1)-(H4) from the previous section hold and we use the same notation for the time and space discretizations. For a two-point approximation of the discrete gradients, we require additionally that the mesh is admissible in the sense of [12, Definition 9.1]. This means that a family of points $(x_{K})_{K\in{\mathcal{T}}}$ is associated to the cells and that the line connecting the points $x_{K}$ and $x_{L}$ of two neighboring cells $K$ and $L$ is perpendicular to the edge $K|L$. For $\sigma\in{\mathcal{E}}_{\rm int}$ with $\sigma=K|L$, we denote by ${\mathrm{d}}_{\sigma}={\mathrm{d}}(x_{K},x_{L})$ the Euclidean distance between $x_{K}$ and $x_{L}$, while for $\sigma\in{\mathcal{E}}_{\rm ext}$, we set ${\mathrm{d}}_{\sigma}={\mathrm{d}}(x_{K},\sigma)$. For a given edge $\sigma\in{\mathcal{E}}$, the transmissibility coefficient is defined by $\tau_{\sigma}=\text{m}(\sigma)/{\mathrm{d}}_{\sigma}$, where $\text{m}(\sigma)$ denotes the Lebesgue measure of $\sigma$.

For the definition of the scheme, we approximate the initial, boundary, and given functions on the elements $K\in{\mathcal{T}}$ and edges $\sigma\in{\mathcal{E}}$:

[TABLE]

and we set $u_{0,K}^{I}=1-\sum_{i=1}^{n}u_{i,K}^{I}$ and $\overline{u}_{0,\sigma}=1-\sum_{i=1}^{n}\overline{u}_{i,\sigma}$. Furthermore, we introduce the discrete gradients

[TABLE]

The numerical scheme is now defined as follows. Let $K\in{\mathcal{T}}$, $k\in\{1,\ldots,N\}$, $i\in\{1,\ldots,n\}$, and $u_{i,K}^{k-1}\geq 0$ be given. Then the values $u_{i,K}^{k}$ are determined by the implicit Euler scheme

[TABLE]

where the fluxes ${\mathcal{F}}_{i,K,\sigma}^{k}$ are given by the upwind scheme

[TABLE]

Here, we have set

[TABLE]

and ${\mathcal{V}}_{i,K,\sigma}^{k}$ is the “drift part” of the flux,

[TABLE]

for $i=1,\ldots,n$. Observe that we employed a double upwinding: one related to the electric potential, defining $\widehat{u}_{0,\sigma,i}^{k}$, and another one related to the drift part of the flux, ${\mathcal{V}}_{i,K,\sigma}^{k}$. The potential is computed via

[TABLE]

The finite volume scheme preserves the structure of the continuous equations only under certain assumptions: {labeling}(A44)

$\partial\Omega=\Gamma_{N}$, i.e., we impose no-flux boundary conditions on the whole boundary.

The diffusion constants are equal, $D_{i}=D>0$ for $i=1,\ldots,n$.

The drift terms are set to zero, $\Phi\equiv 0$.

Without these assumptions, we can only assure the nonnegativity of the discrete concentrations $u_{i}$, $i=1,\ldots,n$. Since we lack a maximum principle for cross-diffusion systems, the upper bounds can only be proven if we assume equal diffusion constants (A2). Under this assumption, the solvent concentration satisfies

[TABLE]

for which a discrete maximum principle can be applied. The $L^{\infty}$ bounds on the concentrations then ensure the existence of solutions for the scheme. If additionally the drift term vanishes (A3), a discrete version of the entropy inequality, the uniqueness of discrete solutions and most importantly, the convergence of the scheme can be proven (under an additional regularity assumption on the mesh). For details, we refer to [5].

4. Numerical experiments

4.1. Implementation

The finite-element discretization is implemented within the finite-element library NGSolve/Netgen, see [30, 31]. The nonlinear equations are solved in every time step by Newton’s method in the variables $w_{i}$ and $\Phi$. The Jacobi matrix is computed using the NGSolve function AssembleLinearization. The finite-volume scheme is implemented in Matlab. Also here, the nonlinear equations are solved by Newton’s method in every time step, using the variables $u$, $\Phi$, and $u_{0}$.

We remark that the finite-volume scheme also performs well when we use a simpler semi-implicit scheme, where we compute $u$ from equation (27) with $\Phi$ taken from the previous time step via Newton’s method and subsequently only need to solve a linear equation to compute the update for the potential. It turned out that this approach is not working for the finite-element discretization. Furthermore, the computationally cheaper implementation used in [23] for a similar scheme in one space dimension, where a Newton and Picard iteration are combined, did not work well in the two-dimensional test cases presented in this paper.

4.2. Test case 1: calcium-selective ion channel

Our first test case models the basic features of an L-type calcium channel (the letter L stands for “long-lasting”, referring to the length of activation). This type of channel is of great biological importance, as it is present in the cardiac muscle and responsible for coordinating the contractions of the heart [6]. The selectivity for calcium in this channel protein is caused by the so-called EEEE-locus made up of four glutamate residues. We follow the modeling approach of [28], where the glutamate side chains are each treated as two half charged oxygen ions, accounting for a total of eight $O^{1/2-}$ ions confined to the channel. In contrast to [28], where the oxygen ions are described by hard spheres that are free to move inside the channel region, we make a further reduction and simply consider a constant density of oxygen in the channel that decreases linearly to zero in the baths (see Figure 1),

[TABLE]

where the scaled maximal oxygen concentration equals $u_{{\rm ox},\max}=(N_{A}/u_{\rm typ})\cdot 52\,\,$mol/L. Here, $N_{A}\approx 6.022\cdot 10^{23}\,$mol*-1* is the Avogadro constant and $u_{\rm typ}=3.7037\cdot 10^{25}L^{-1}$ the typical concentration (taken from [4, Table 1]). In addition to the immobile oxygen ions, we consider three different species of ions, whose concentrations evolve according to model equations (1): calcium (Ca2+, $u_{1}$), sodium (Na+, $u_{2}$), and chloride (Cl-, $u_{3}$). We assume that the oxygen ions not only contribute to the permanent charge density $f=-u_{\rm ox}/2$, but also take up space in the channel, so that we have $u_{0}=1-\sum_{i=1}^{3}u_{i}-u_{\rm ox}$ for the solvent concentration.

For the simulation domain, we take a simple geometric setup resembling the form of a channel; see Figure 1. The boundary conditions are as described in the introduction, with constant values for the ion concentrations and the electric potential in the baths. The physical parameters used in our simulations are taken from [4, Table 1]. The simulations are performed with a constant (scaled) time step size $\tau=2\cdot 10^{-4}$. The initial concentrations are simply taken as linear functions connecting the boundary values. An admissible mesh consisting of 74 triangles was created with Matlab’s initmesh command, which produces Delauney triangulations. Four finer meshes were obtained by regular refinement, dividing each triangle into four triangles of the same shape.

We remark that the same test case was already used in [5] to illustrate the efficiency of the finite-volume approximation. Furthermore, numerical simulations for a one-dimensional approximation of the calcium channel can be found in [4] for stationary solutions and in [15] for transient solutions.

Figures 2 and 3 present the solution to the ion transport model in the original variables $u$ and $\Phi$ at two different times; the first one after only 600 time steps and the second one after 6000 time steps, which is already close to the equilibrium state. The results are computed on the finest mesh with 18,944 elements. In the upper panel, the concentration profiles and electric potential as computed with the finite-element scheme are depicted. In the lower panel, the difference between the finite-volume and finite-element solutions is plotted. We have omitted the plots for the third ion species (Cl-), since it vanishes almost immediately from the channel due to its negative charge. While absolute differences are relatively small, we can still observe that the electric potential in the finite-element case is always smaller compared to the finite-volume solution, while the peaks of the concentration profiles are more distinctive for the finite-element than for the finite-volume solution.

In order to compare the two numerical methods, we test the convergence of the schemes with respect to the mesh diameter. Since an exact solution to our problem is not available, we compute a reference solution both with the finite-volume and the finite-element scheme on a very fine mesh with 18,944 elements and maximal cell diameter $h\approx 0.01$. The differences between these reference solutions in the discrete $L^{1}$ and $L^{\infty}$ norms are given in Table 1 for the various unknowns. Since the finite-element and finite-volume solutions are found in different function spaces, one has to be careful how to compare them. The values in Table 1 are obtained by projecting the finite-element solution onto the finite-volume space of functions that are constant on each cell in NGSolve, thereby introducing an additional error. However, the difference between the reference solutions is still reasonably small, especially when the simulations are already close to the equilibrium state.

To avoid the interpolation error in the convergence plots, we compare the approximate finite-element or finite-volume solutions on coarser nested meshes with the reference solutions computed with the corresponding method. In Figure 4, the errors in the discrete $L^{1}$ norm between the reference solution and the solutions on the coarser meshes at the two fixed time steps $k=600$ and $k=6000$ are plotted. For the finite-volume approximation, we clearly observe the expected first-order convergence in space, whereas for the finite-element method, the error decreases, again as expected, with $h^{2}$. These results serve as a validation for the theoretical convergence result proven for the finite-element scheme and show the efficiency of the finite-volume method even in the general case of ion transport, which is not covered by the convergence theorem in [5].

In Table 2, the average time needed to compute one time step with the finite-element or finite-volume scheme for the five nested meshes is given. Clearly, the finite-volume scheme is much faster than the finite-element method. This is mostly due to the computationally expensive assembly of the finite-element matrices.

In addition to the convergence analysis, we also study the behavior of the discrete entropy for both schemes. We consider in both cases the entropy relative to the steady state $(u^{\infty}_{i},\Phi^{\infty})$, which is computed from the corresponding discretizations of the stationary equations with the same parameters and boundary data. Figure 5 shows the relative entropy (see [5, Section 6]) and the $L^{1}$ error compared to the equilibrium state for the finite-element and finite-volume solutions on different meshes. Whereas for the coarsest mesh the convergence rates differ notably, we can observe a similar behavior when the mesh is reasonably fine. In Figure 6, we investigate the convergence of the relative entropy with respect to the mesh size. As before, we observe second-order convergence for the finite-element scheme and a first-order rate for the finite-volume method.

4.3. Test case 2: bipolar ion channel

The second example models a pore with asymmetric charge distribution, which occurs naturally in biological ion channels but also in synthetic nanopores. Asymmetric pores typically rectify the ion current, meaning that the current measured for applied voltages of positive sign is higher than the current for the same absolute value of voltage with negative sign. The setup is similar to that of an N-P semiconductor diode. The N-region is characterized by the fixed positive charge. The anions are the counter-ions and thus the majority charge carriers, while the cations are the co-ions and minority charge carriers. In the P-region, the situation is exactly the other way around. In the on-state, the current is conducted by the majority carriers, while in the off-state, the minority carriers are responsible for the current, which leads to the rectification behavior.

Often, bipolar ion channels are modeled with asymmetric surface charge distributions on the channel walls. However, to fit these channels into the framework of our model, we follow the approach described in [17]. Similar to the first test case, we assume that there are eight confined molecules inside the channel, but this time four molecules are positively charged ($+0.5e$) and the other four molecules are negatively charged ($-0.5e$). The simulation domain $\Omega\subset{\mathbb{R}}^{2}$ is depicted in Figure 7. The shape of the domain and the parameters used for the simulations are taken from [17] and are summarized in Table 3. The mesh (made up of 2080 triangles) was created with NGSolve/Netgen. We consider two mobile species of ions, one cation (Na+, $u_{1}$) and one anion (Cl-, $u_{2}$). The confined ions are modeled as eight fixed circles of radius $1.4$, where the concentration $c\equiv c_{\max}$ is such that the portion of the channel occupied by these ions is the same as in the simulations in [17]. The solvent concentration then becomes $u_{0}=1-u_{1}-u_{2}-c$.

By changing the boundary value $\overline{\Phi}_{\text{right}}$ for the potential $\Phi$ on the right part of the Dirichlet boundary (on the left side, it is fixed to zero), we can apply an electric field in forward bias (on-state, $\overline{\Phi}_{\text{right}}=1$) or reverse bias (off-state, $\overline{\Phi}_{\text{right}}=-1$). Figures 8 and 9 show the stationary state computed with the finite-element method in the on- and off-state, respectively. Evidently, the ion concentrations in the on-state are much higher than in the off-state. In comparison with the results from [17], where the Poisson–Nernst–Planck equations with linear diffusion (referred to as the linear PNP model) were combined with Local Equilibrium Monte–Carlo simulations, we find that with the Poisson–Nernst–Planck equations with cross-diffusion (referred to as the nonlinear PNP model), the charged ions in the channel attract an amount of ions higher than the bath concentrations even in the off-state.

From a modeling point of view, it is an important question whether the nonlinear PNP model reproduces the rectification mechanism described above. For this purpose, we need to calculate the electric current $I$ flowing through the pore, given by

[TABLE]

where $A$ is the cross-section of the pore and $\nu$ the unit normal to $A$. In the finite-element setting, we can use the representation of the fluxes in entropy variables, ${\mathcal{F}}_{i}=D_{i}u_{i}(w,\Phi)u_{0}(w,\Phi)\nabla w_{i}$ and compute the integrals in (28) using a quadrature formula along the line $x=10$.

Figure 10 shows the current-voltage curves obtained with the finite-element solutions. In addition, the rectification is depicted, which is calculated for voltages $U\geq 0$ according to

[TABLE]

We also compute the current-voltage curve for the linear PNP model, which is obtained from the model equations by setting $u_{0}\equiv 1$, such that

[TABLE]

We expect from the simulations done in [4] for the calcium channel that the current of the nonlinear PNP model is lower than that one from the linear PNP model. This expectation is confirmed also in this case. As Figure 10 shows, the rectification is stronger in the nonlinear PNP model. The difference between the two models is even more pronounced when we increase the concentration of the confined ions to $c_{\max}=0.7$. In that case, the channel gets more crowded and size exclusion has a bigger effect. We observe a significantly lower current and higher rectification for the nonlinear PNP model.

5. Conclusions

In this work, we have presented a finite-element discretization of a cross-diffusion Poisson–Nernst–Planck system and recalled a finite-volume scheme that was previously proposed and analyzed [5]. In the following, we summarize the differences between both approaches from a theoretical viewpoint and our findings from the numerical experiments.

• •

Structure of the scheme: The finite-element scheme strongly relies on the entropy structure of the system and is formulated in the entropy variables. From a thermodynamic viewpoint, the entropy variables are related to the chemical potentials, which gives a clear connection to nonequilibrium thermodynamics. On the other hand, the finite-volume scheme exploits the drift-diffusion structure that the system displays in the original variables.

• •

$L^{\infty}$** bounds:** Due to the formulation in entropy variables, the $L^{\infty}$ bounds for the finite-element solutions follow immediately from (8) without the use of a maximum principle. In other words, the lower and upper bounds are inherent in the entropy formulation. In the case of the finite-volume scheme, we can apply a discrete maximum principle, but only under the (restrictive) assumption that the diffusion coefficients $D_{i}$ are the same.

• •

Convergence analysis: The entropy structure used in the finite-element scheme allows us to use the same mathematical techniques for the convergence proof as for the continuous system, but a regularizing term has to be added to ensure the existence of discrete solutions. The convergence of the finite-volume solution requires more restrictive assumptions: In addition to the equal diffusion constants necessary for proving the existence and $L^{\infty}$ bounds, we can only obtain the entropy inequality and gradient estimates for vanishing potentials.

• •

Initial data: Since the initial concentrations have to be transformed to entropy variables via (6), the finite-element scheme can only be applied for initial data strictly greater than zero. The finite-volume scheme, on the other hand, can handle exactly vanishing initial concentrations.

• •

Experimental convergence rate: In the numerical experiments, both schemes exhibit the expected order of convergence with respect to mesh size (even if we cannot prove any rates analytically): first-order convergence for the finite-volume scheme and second-order convergence for the finite-element scheme.

• •

Performance: The numerical experiments done for this work suggest that the finite-element algorithm needs smaller time steps for the Newton iterations to converge than for the finite-volume scheme, especially when the solvent concentration is close to zero. Furthermore, the assembly of the finite-element matrices is computationally quite expensive resulting in longer running times compared to the finite-volume scheme.

• •

Mesh requirements: A finite-volume mesh needs to satisfy the admissibility condition. This might be a disadvantage for simulations in three space dimensions.