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

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.


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 (Hodgkin and Huxley 1952) as well as Brownian and Langevin dynamics (Im et al. 2000;Nadler et al. 2005) to the widely used Poisson-Nernst-Planck model (Eisenberg 1998). 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 Communicated by Jorge X. Velasco. B Ansgar Jüngel juengel@tuwien.ac.at Anita Gerstenmayer anita.gerstenmayer@tuwien.ac.at terms (Gillespie et al. 2005;Nonner et al. 2000), 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 (Burger et al. 2010;Simpson et al. 2009). These nonlinear cross-diffusion terms are common in diffusive multicomponent systems (Jüngel 2016, Chapter 4). A lattice-free approach, starting from stochastic Langevin equations, can be found in Bruna and Chapman (2014). The scope of this paper is to present a new finiteelement discretization of the degenerate cross-diffusion system and to compare this scheme to a previously proposed finite-volume method (Cancès et al. 2019).
The dynamics of the ion concentrations u = (u 1 , . . . , u n ) is governed by the evolution equations where u 0 = 1 − n i=1 u i denotes the solvent concentration, D i > 0 is the diffusion constant, z i the ion charge, and β a mobility parameter. To be precise, u i is the mass fraction of the ith ion, and we refer to n i=0 u i = 1 as the total relative mass, just meaning that the ionsolvent mixture is saturated. The electric potential is self-consistently given by the Poisson equation with the permanent charge density f = f (x) and the scaled permittivity constant λ 2 . The equations are solved in a bounded domain ⊂ R d with smooth boundary ∂ . Equations (1) are equipped with initial data u(0) = u I satisfying 0 < n i=1 u I i < 1. The boundary ∂ consists of an insulating part N and the union D of contacts with external reservoirs: System (1)-(2) can be interpreted as a generalized Poisson-Nernst-Planck model. The usual Poisson-Nernst-Planck equations (Eisenberg 1998) follow from (1) by setting u 0 = 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) (Hsieh et al. 2015). 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 Dreyer et al. (2013), but the coefficients differ from (1). The model at hand was derived in Burger et al. (2010) and Simpson et al. (2009) from a lattice model taking into account finite-size effects.
Model (1)-(4) contains some mathematical difficulties. First, its diffusion matrix A(u) = (A i j (u)) ∈ R n×n , given by A i j (u) = D i u i for i = j and A ii (u) = D i (u 0 + u i ) for i = 1, . . . , 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 ≤ u i ≤ 1 for i = 1, . . . , 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, Indeed, writing Eq. (1) in terms of the entropy variables w 1 , . . . , w n , given by it follows that where the new diffusion matrix is symmetric and positive semidefinite (in fact, it is even diagonal). This procedure has a thermodynamical background: the quantities ∂h/∂u i are known as the chemical potentials, and B is the so-called mobility or Onsager matrix [de Groot and Mazur (1984)]. The transformation to entropy variables also solves the third difficulty. Solving the transformed system (7) for w = (w 1 , . . . , w n ), the concentrations are given by showing that u i is positive and bounded from above: Moreover, the entropy structure leads to gradient estimates via the entropy inequality 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 using the "degenerate" Aubin-Lions lemma of Jüngel 2015, Appendix C or its discrete version in Cancès et al. 2019, Lemma 10. These ideas were employed in Burger et al. (2010) 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 Jüngel (2015), Zamponi and Jüngel (2017), still excluding the potential. A global existence result for the full problem (1)-(4) was established in Gerstenmayer and Jüngel (2018).
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 Cancès et al. (2019), analyzing a finitevolume 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 Lu et al. (2010), Prohl and Schmuck (2009), there seems to be no numerical analysis of the ion-transport model (1)-(4) apart from the finite-volume scheme in Cancès et al. (2019) and simulations of the stationary equations in Burger et al. (2012). Let us mention some other works on finite-element methods for related cross-diffusion models. In Barrett and Blowey (2004), 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 also employed in Galiano and Selgas (2014). A lumped finite-element method was analyzed in Frittelli et al. (2017) for a reaction-cross-diffusion equation on a stationary surface with positive definite diffusion matrix. In Jüngel and Leingang (2018), 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 Egger (2018). 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 (Sect. 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 the 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 Cancès et al. (2019) (recalled in Sect. 3) are applied to two test cases in two space dimensions: a calcium-selective ion channel and a bipolar ion channel (Sect. 4). Static current-voltage curves show the rectifying behavior of the bipolar ion channel. • The advantages and drawbacks of both schemes are discussed (Sect. 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.

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: (H1) Domain: (H4) Initial and boundary data: 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 τ > 0 and set t k = kτ for k = 1, . . . , N , where N ∈ N is given and τ = T /N .
For the space discretization, we introduce a family T h (h > 0) of triangulations of , consisting of open polygonal convex subsets of (the so-called cells) such that = ∪ K ∈T h K with maximal diameter h = max K ∈T h diam(K ). We assume that the corresponding family of edges E can be split into internal and external edges, Each exterior edge is assumed to be an element of either the Dirichlet or Neumann boundary, i.e., E ext = E D ext ∪ E N ext . For given K ∈ T h , we define the set 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 In the finite-element setting, the triangulation is completed by the set of nodes { p j : j ∈ J }. We impose the following regularity assumption on the mesh. There exists a constant γ ≥ 1 such that where ρ K is the radius of the incircle and h K is the diameter of K .
We associate with T h the usual conforming finite-element spaces and H 1 D ( ) is the set of H 1 ( ) functions that vanish on D in the weak sense. Let {χ j } j∈J be the standard basis functions for S(T h ) with χ j ( p i ) = δ i j for all i, j ∈ J . We define the nodal interpolation operator I h : Due to the regularity assumptions on the mesh, I h has the following approximation property [(see, e.g., (Ciarlet 1978 (10)

Definition of the scheme
To define the finite-element scheme, we need to approximate the initial and boundary data. We set w 0 where 0 is the standard finite-element solution to the linear equation (2) with u I i on the right-hand side. Furthermore, we set The finite-element scheme is now defined as follows. Given for all φ ∈ S D (T h ) n and θ ∈ S D (T h ). The symbol ":" signifies the Frobenius matrix product; here, the expression reduces to The term involving the parameter ε > 0 is only needed to guarantee the coercivity of (11), (12). Indeed, the diffusion matrix B(w k , k ) degenerates when w k i → −∞, and the corresponding bilinear form is only positive semidefinite. To emphasize the dependence on the mesh and ε, we should rather write w (h,ε,k) instead of w k and similarly for k ; however, for the sake of presentation, we will mostly omit the additional superscripts. The original variables are recovered by computing . . . , N , and u (τ ) (·, 0) = I h u I as well as similarly for (τ ) , we obtain piecewise constant in time functions.

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: where H is defined in (5) and The proof of the lemma is similar to the proof of Theorem 1 in Gerstenmayer and Jüngel (2018). The main difference is that in Gerstenmayer and Jüngel (2018), a regularization term of the type ε (− ) m w k + w k has been added to achieve via H m ( ) → L ∞ ( ) for m > d/2 compactness and L ∞ solutions. In the finite-dimensional setting, this embedding is not necessary but we still need the regularization ε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 ε in the numerical simulations such that the additional term does not affect the solution practically.

Proof of Lemma 1
Let y ∈ S(T h ) n and δ ∈ [0, 1]. There exists a unique solution k to (12) with w k replaced by y+w h , satisfying k − h ∈ S D (T h ), since the function → z i u i (y, ) is bounded and nonincreasing. Indeed, let 1 and 2 be two solutions with the same initial datum. Then, taking the difference of the corresponding weak formulations (12) and using the test function 1 − 2 , we find that holds for some constant C > 0. Next, we consider the linear problem where The bilinear form a and the linear form F are continuous on S D (T h ) n . The equivalence of all norms on the finite-dimensional space S D (T h ) implies the coercivity of a, shows that all elements v are bounded independent of y and δ and thus, all fixed points v = S(v, δ) are uniformly bounded. Furthermore, S(y, 0) = 0 for all y ∈ S D (T h ) n . The continuity of S follows from standard arguments and the compactness comes from the fact The discrete entropy inequality (13) is proven using τ (w k − w h ) ∈ S D (T h ) n as a test function in (11) and exploiting the convexity of H , which concludes the proof.
Remark 1 (Structure preservation of the scheme) Lemma 1 shows that if the boundary data are in thermal equilibrium, i.e., ∇w h = 0, then the finite-element scheme (11), (12) dissipates the entropy (5), i.e., H u k ≤ H u k−1 . Moreover, it preserves the invariant region D, i.e., u k ∈ D, and the mass fraction u k i is nonnegative and bounded by one. The scheme conserves the total relative mass, i.e., n i=1 u k i L 1 ( ) = 1, which is a direct consequence of the definition of u k 0 .

Uniform estimates
The next step is the derivation of a priori estimates uniform in the parameters ε, τ , 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: for i = 1, . . . , n, where here and in the following, C > 0 is a generic constant independent of ε, τ , and h.
Proof As the proof is similar to that one in the continuous setting, we give only a sketch. Note that our finite-element space is a subset of H 1 ( ), so the computations can be done as in Gerstenmayer and Jüngel (2018), and in particular, the chain rule holds. Observe that the definition of the entropy variables implies that 0 < u k i < 1 in for i = 1, . . . , n and k = 1, . . . , N . It is shown in the proof of Lemma 6 of Gerstenmayer and Jüngel (2018) that (13) gives We resolve this recursion to find that The right-hand side is bounded because of (14), τ k ≤ T , and the boundedness of the interpolation operator. Inserting the identity the estimates follow.

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).
Proof We pass first to the limit (ε, h) → 0 and then τ → 0, since the latter limit can be performed as in the proof of Theorem 1 in Gerstenmayer and Jüngel (2018). Fix k ∈ {1, . . . , N } and let u ε,h,k) be the approximate solution from Lemma 1. We set u . Using the compact embedding H 1 ( ) → L 2 ( ) and the a priori estimates from Lemma 2, it follows that there exists a subsequence which is not relabeled such that, as (ε, h) → 0, Combining (17) and (21), we infer that (up to a subsequence) Next, let φ ∈ (H 2 ( ) ∩ H 1 D ( )) n . As we cannot use φ i directly as a test function in (11), we take I h φ ∈ S D (T h ) n , where I h is the interpolation operator, see (10). To pass to the limit in (11), we rewrite the integral involving the diffusion matrix: (ε,h) )∇w (ε,h) We estimate each of the above summands separately. For the last term, we proceed as follows: Similarly as for (24), it follows that Then, together with the weak convergence of ∇ (ε,h) , we infer that h) ) is bounded in L 2 ( ), this weak convergence also holds in L 2 ( ). Because of the interpolation property (10) and estimate (26), Following the arguments of Step 3 in Gerstenmayer and Jüngel (2018), Section 2, [using (24)], we have Thus, the limit (ε, h) → 0 in (25) gives (ε,h) , (ε,h) )∇w (ε,h) Furthermore, we deduce from (23) that Thus, passing to the limit (ε, h) → 0 in scheme (11), (12) leads to A density argument shows that we can take test functions φ i , θ ∈ H 1 D ( ). The a priori estimates from Lemma 2 remain valid in the weak limit. Now the limit τ → 0 can be done exactly as in Gerstenmayer and Jüngel (2018), Theorem 1, Step 4, which concludes the proof.

The finite-volume scheme
We briefly recall the finite-volume scheme from Cancès et al. (2019) and summarize the assumptions and results, as this is necessary for the comparison of the finite-element and finite-volume scheme in Sect. 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 Eymard et al. (2000), Definition 9.1. This means that a family of points (x K ) K ∈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 σ ∈ E int with σ = K |L, we denote by d σ = d(x K , x L ) the Euclidean distance between x K and x L , while for σ ∈ E ext , we set d σ = d(x K , σ ). For a given edge σ ∈ E, the transmissibility coefficient is defined by τ σ = m(σ )/d σ , where m(σ ) denotes the Lebesgue measure of σ .
For the definition of the scheme, we approximate the initial, boundary, and given functions on the elements K ∈ T and edges σ ∈ E: and we set u I 0,K = 1 − n i=1 u I i,K and u 0,σ = 1 − n i=1 u i,σ . Furthermore, we introduce the discrete gradients The numerical scheme is now defined as follows. Let K ∈ T , 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 F k i,K ,σ are given by the upwind scheme Here, we have set and V 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 u k 0,σ,i , and another one related to the drift part of the flux, V k i,K ,σ . The potential is computed via The finite-volume scheme preserves the structure of the continuous equations only under certain assumptions: Without these assumptions, we can only assure the nonnegativity of the discrete concentrations u i , i = 1, . . . , 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 for which a discrete maximum principle can be applied. The L ∞ 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 Cancès et al. (2019).

Implementation
The finite-element discretization is implemented within the finite-element library NGSolve/ Netgen, see Schöberl (1997Schöberl ( , 2014. The nonlinear equations are solved in every time step by Newton's method in the variables w i and . 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, , and u 0 . The integrals appearing in scheme (11), (12) are computed using a Gauß quadrature implemented in NGSolve that computes the trial functions exactly. Because of the nonlinear functions appearing in the integrals, the quadrature yields only approximate values.
We remark that the finite-volume scheme also performs well when we use a simpler semiimplicit scheme, where we compute u from Eq. (27) with 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 Jüngel and Leingang (2018) 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.

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 (Carafoli et al. 2001). 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 Nonner et al. (2001), 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 Nonner et al. (2001), 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 Fig. 1 where the scaled maximal oxygen concentration equals u ox,max = (N A /u typ ) × 52 mol/L. Here, N A ≈ 6.022 × 10 23 mol −1 is the Avogadro constant and u typ = 3.7037 × 10 25 L −1 the typical concentration [taken from (Burger et al. 2012, 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 (Ca 2+ , 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 ox /2, but also take up space in the channel, so that we have u 0 = 1 − 3 i=1 u i − u ox for the solvent concentration.
For the simulation domain, we take a simple geometric setup resembling the form of a channel; see Fig. 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 (Burger et al. 2012, Table 1). The simulations are performed with a constant (scaled) time step size τ = 2 × 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 Cancès et al. (2019) 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 Burger et al. (2012) for stationary solutions and in Gerstenmayer and Jüngel (2018) for transient solutions. Figures 2 and 3 present the solution to the ion-transport model in the original variables u and 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.
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 ≈ 0.01. The differences between these reference solutions in the discrete L 1 and L ∞ 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 Fig. 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 Cancès et al. (2019).
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 ∞ i , ∞ ), 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 (Cancès et al. 2019, 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 Fig. 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.

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 Ható et al. (2016). 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 ⊂ R 2 is depicted in Fig. 7. The shape of the domain and the parameters used for the simulations are taken from Ható et al. (2016) 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 ≡ c max is such that the portion of the channel occupied by these ions is the same as in the simulations in Ható et al. (2016). The solvent concentration then becomes u 0 = 1 − u 1 − u 2 − c.
By changing the boundary value right for the potential 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 (onstate, right = 1) or reverse bias (off-state, 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 Ható et al. (2016), 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 crossdiffusion (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 where A is the cross section of the pore and ν the unit normal to A. In the finiteelement setting, we can use the representation of the fluxes in entropy variables, F i = D i u i (w, )u 0 (w, )∇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 ≥ 0 according to .
We also compute the current-voltage curve for the linear PNP model, which is obtained from the model equations by setting u 0 ≡ 1, such that ∂ t u i = div D i ∇u i + D i βz i u i ∇ .
We expect from the simulations done in Burger et al. (2012) 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 also confirmed in this case. As Fig. 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  Table 3; second row: with c max = 0.7 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.

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 (Cancès et al. 2019). 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 ∞ bounds Due to the formulation in entropy variables, the L ∞ bounds for the finiteelement 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 ∞ 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 finiteelement 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.