Numerical analysis of a finite element formulation of the P2D model for Lithium-ion cells

The mathematical P2D model is a system of strongly coupled nonlinear parabolic-elliptic equations that describes the electrodynamics of lithium-ion batteries. In this paper, we present the numerical analysis of a finite element-implicit Euler scheme for such a model. We obtain error estimates for both the spatially semidiscrete and the fully discrete systems of equations, and establish the existence and uniqueness of the fully discrete solution.


Introduction
In this paper, we present the numerical analysis of a finite element-implicit Euler method to calculate the numerical solution of the so called pseudo-two-dimensional (P2D) model. proposed by J. Newman and coworkers [3]. This is a mathematical model based on the electrochemical kinetics and continuun mechanics laws, which consists of a system of coupled nonlinear parabolic-elliptic equations to model the physical-chemical phenomena governing the behavior of lithium ion batteries. The P2D model is very much used in engineering studies. A good presentation of it can be found in [13] and [15]. A lithium-ion battery system is composed of a number of lithium-ion cells. A typical cell consists of three regions, namely, a porous negative electrode (which plays the role of anode of the cell in the discharge process) connected to the negative terminal collector of the battery, a separator that is an electron insulator allowing the flow of lithium ions between the anode and the cathode, and a porous positive electrode (which plays the role of cathode during the discharge process) connected to the positive terminal, see Fig. 1. We must point out that in the charge process the negative electrode plays the role of cathode and the positive electrode is the anode. The electrodes are composite porous structures of highly packed active lithium particles, typically Li x C 6 in the negative electrode and metal oxide, such as Li 1−x Mn 2 O 4 , in the positive electrode, plus a binder and a polymer that act as conductive agents. Furthermore, the cell is filled with the electrolyte that occupies the holes left free by the particles and the filler material. The electrolyte is a lithium salt dissolved in an organic solvent. In the description of the model it is customary to consider two phases: the electrolyte phase and the solid phase, the latter is composed of the solid particles of the electrodes.
The P2D model of a lithium-ion cell considers that the dynamics is only relevant along the x-axis, neglecting what happens along the y-axis and z-axis, because the ratios Lx Ly and Lx Lz = O(10 −3 ), L x , L y and L z being the characteristic length scales along the corresponding axes. The main modeling assumptions are the following : (1) The active particles of the electrodes are assumed to be spheres of radius R s which may be different in each electrode. (2) Side reactions are neglected and no gas phase is present. (3) The transport of lithium ions is due to diffusion and migration in the electrolyte solution, and in the solid particles the atoms of lithium move between vacancies in the crystalline structure of the particles due to local diffusion in concentration. By longitudinal and latitudinal symmetry considerations, the diffusion in the active particles is only in the radial direction. (4) The electrochemical reaction of lithium insertion and extraction processes follows the Buttler-Volmer law. (5) The effective transport coefficients are calculated by the Bruggeman relation, i.e., µ ef f = µε p (p=1.5), where µ is a generic where D n , D s and D p denote the domains of the negative electrode, the separator and the positive electrode respectively. Notice that D 1 represents the cell domain, D 2 is a domain that is the union of two disjoint domains corresponding to the electrodes, and D 3 is in a certain sense a modeling space accounting for the spherical balls of radius R s (x) that represent at each x ∈ D 2 the solid active particles, such that when x ∈ D n , R s (x) = R − s , and when x ∈ D p , R s (x) = R + s . The variables of the model are the following: for the electrolyte phase, the molar concentration of lithium ions u(x, t), and the electric potential φ 1 (x, t), x ∈ D 1 ; for the solid phase, the molar concentration of lithium v(x; r, t), x ∈ D 2 and r ∈ (0; R s (x)), and the electric potential φ 2 (x, t), x ∈ D 2 . Another important variable is the so called molar flux of lithium ions exiting the solid particles, J(x, u, v, φ 1 , φ 2 , U )/F , F being the Faraday constant. The mathematical expression of J is given by the Buttler-Volmer law, see (1).
Many numerical models to integrate the P2D model have been proposed. The first one is the Dualfoil model developed by J. Newman and his collaborators [12], this is a model that uses second order finite differences for space discretization of the differential operators combined with the first order backward Euler time stepping scheme; the Dualfoil model is distributed as free software, which is being updated through time. Later on, authors such as [11] and [17], just to cite a few, have developed their own codes by using second order finite volume for space discretizations combined with the first order in time implicit Euler scheme for time discretization. Other authors make the numerical simulations with COMSOL multi-physics package that uses finite elements for space discretizations of the equations, the resulting system of nonlinear differential equations is integrated by different time stepping schemes, in particular, conventional DAE solvers, such as DASK [14]. New numerical models have recently been proposed to improve the computational efficiency, to this respect, we mention the operator splitting technique of [6], the orthogonal collocation method for space discretization combined with the first order implicit Euler scheme for time discretization of [9], and the implicit-explicit Runge-Kutta-Chebyshev finite element method of [1]. Despite the activity in the development of numerical methods no rigorous numerical analysis of such methods has been published so far; so, to the best of our knowledge, this is the first paper presenting the analysis of a numerical method developed to integrate the P2D model.
The layout of the paper is the following. In Section 2 we introduce the governing equations of the P2D model together with the functional framework needed for the numerical analysis. Section 3 is devoted to the semidiscrete space discretization of the model in a finite element framework. The error analysis of the semi-discrete solution is performed in Section 4. Since this analysis is long, then we have split the section into three subsections in order to make more palatable its presentation. Subsection 4.1 is a collection of auxiliary results; subsections 4.2 and 4.3 deal with the error estimates for the potentials and the concentrations, respectively. The fully discrete model and its error analysis is presented in Section 5, which is also split into subsections. Since the fully discrete model is a nonlinear system of elliptic and fully discrete parabolic equations at each time instant t n , then we have also studied the existence and uniqueness of the solution by applying Minty-Browder theorem [18] for the elliptic equations, and Brower´s fixed point theorem for the parabolic equations.

The governing equations of the isothermal P2D model
We consider the governing equations of the isothermal P2D model for the variables u(x, t), v(x; r, t), φ 1 (x, t) and φ 2 (x, t) presented in Chapters 3 and 4 of [15]. However, to facilitate both the formulation of the numerical method to integrate these equations and its numerical analysis, it is convenient to make the changes of variable introduced in [10] and [19]. Thus, in order to make homogeneous the Neumann type boundary conditions for the potential φ 2 one considers the function H(x, t) given by the expression where I(t) denotes the applied current, A is the area of the plate and σ is a positive coefficient defined below, and replace φ 2 (x, t) by φ 2 (x, t) + H(x, t); likewise, we replace the potential , where κ(u) > 0 denotes the effective electrolyte phase ionic conductivity; t 0 + > 0 is the so called transfer number, which is assumed to be constant; R is the universal gas constant and T denotes the absolute temperature inside the cell, which is assumed to be constant in the isothermal model; this latter change of variable for φ 1 (x, t) simplifies the expression of the equation for the potential of the electrolyte phase written in Chapter 4 of [15], making it more manageable from a computational viewpoint. Another important variable, as we mentioned above, is the reaction current density J. The reaction rate is coupled to phase potentials by the Buttler-Volmer kinetic expression.
In this expression, v s = v(x; R s (x), t) denotes the lithium concentration on the surface of the active is the active area per electrode unit volume; ε s (x) denotes the volume fraction of the active material, ε s (x) = ε − s > 0 for x ∈ D n and ε s (x) = ε + s > 0 for x ∈ D p ; α a ∈ (0, 1) and α c ∈ (0, 1) are anodic and cathodic transfer coefficients for an electron reaction; R SEI represents the solid interface resistance, usually, R SEI = 0 in the engineering literature unless the model also considers aging phenomena of the battery, so in this paper we take R SEI = 0.
where U stands for the equilibrium potential at the solid electrolyte interface, which is assumed to be known. i 0 is the exchange current density, i.e., here, v max is the maximum concentration of lithium in the solid phase, which may have different values in the positive and negative electrodes, so Considering the above mentioned changes of variable and taking the transfer coefficients α a and α c equal to 0.5, as many engineering papers do, the expression for the reaction current that we use in the paper is where β = F 2RT ; a 2 (x) = 3ε s (x)/R s (x); and Noting that the boundaries ∂D 1 and ∂D 2 of the domains D 1 and D 2 are ∂D 1 := {0, L} and ∂D 2 := {0, L n , L n + δ, L}, we formulate the equations of the model as follows.
Concentration u(x, t) in the electrolyte phase.
Concentration v(x; r, t) in the solid phase. For almost every Solid phase potential φ 2 (x, t). where . In these equations, k 1 (x) > 0 and k 2 (x) > 0 represent effective diffusion coefficients in the electrolyte and solid phases respectively, and σ(x) denotes the effective electric conductivity in the solid phase. The functions a 1 (x), a 2 (x) , σ(x) and k 2 (x) are considered to be piecewise positive constant functions in the sense that they have different constant values in the negative electrode, separator and positive electrode.
We also have to consider that for t ∈ (0, T end ), J(x, u, v s , φ 1 , φ 2 ) satisfies the algebraic conditions Notice that the first row of algebraic conditions follow directly from (7) and the definition of J, whereas the second row conditions translates the boundary conditions of the solid phase potential. It is worth remarking the conservative properties enjoyed by both u(x, t) and v(x; r, t); namely, for all t ∈ [0, T end ] These relations are readily obtained by integrating (5) and (6) and using the corresponding boundary conditions. Moreover, it can be shown [10] that for t ≥ 0 and x ∈ D 1 , u(x, t) > 0, similarly, for (x, r) ∈ D 3 , 0 < v(x; r, t) < v max . Let D denote a generic open bounded domain in R; hereafter, the closure of a domain D is denoted D. The functional spaces that we use in this paper are the following. The Sobolev spaces H m (D), m being a nonnegative integer, when m = 0, H 0 (D) := L 2 (D); the Lebesgue spaces L p (D), 1 ≤ p ≤ ∞; the spaces of measurable radial functions [4] q being a nonnegative integer, when q = 0 we set H 0 r (0, R) := L 2 r (0, R); also, for p being a nonnegative integer, the normed spaces of measurable functions where v 2 H p (D2;H q r (0,Rs(·))) = p j=0 D2 dx; and the spaces H p,q r (D 2 × (0, R s (·))) = H p (D 2 , L 2 r (0, R s (·))) ∩ L 2 (D 2 , H q r (0, R s (·))) with norm u 2 H p,q r (D2×(0,Rs(·))) := u 2 H p (D2,L 2 r (0,Rs(·))) + u 2 L 2 (D2,H q r (0,Rs(·))) . Notice that H 0,0 r (D 2 × (0, R s (·))) = L 2 (D 2 , L 2 r (0, R s (·))). Since the variables of the model depend on time, then we also introduce the normed spaces L p (0, t; X), where 1 ≤ p ≤ ∞, and (X, · X ) being a real Banach space.
Other spaces used in the paper are W (D 1 ) := v ∈ H 1 (D 1 ) : D1 vdx = 0 , which is a closed subspace of H 1 (D 1 ) where the potential φ 1 (x, t) is calculated, and the space of q times continuously differentiable functions defined on D, C q (D), when q = 0, C 0 (D) := C(D).

Remark 1
Following the arguments of [4], where its non-isothermal P2D model includes an additional time dependent non linear ordinary equation for the bulk temperature T (t), one can formulate an alternative definition of the weak solution to (5)-(8) based on its Definition 2.7 and prove, under the assumptions A1-A3 and for a partition t 0 < t 1 < · · · < t N of [0, T end ], T end being small enough, that there is a unique , v(t n )) being the initial condition in such an interval. Also, Kröner [10] proves a local existence and uniqueness theorem for the weak solution of the isothermal P2D model under less general assumptions than in [4]. 3 The semidiscrete finite element formulation of the isothermal P2D model We use H 1 -conforming linear finite elements (P 1 −finite elements) for the space approximation of the variables u(x, t), φ 1 (x, t) and φ 2 (x, t); however, v(x; r, t) is approximated by nonconforming P 0 −finite elements in the x−coordinate and H 1 -conforming P 1 −finite elements in the r−coordinate. The family of meshes D 1h constructed on the domain D 1 includes the points x = 0, x = L n , x = L n + δ and x = L as mesh points; since these points are also boundary points of D 2 , then they are also considered as mesh points in the family of meshes D 2h . Figure 2 illustrates the families of meshes that we are going to describe next. Noting that D 2 ⊂ D 1 , we choose the family of meshes D 2h as a subset of D 1h . Let N E 1 and N E 2 be the number of elements of D 1h and D 2h respectively, and let M 1 and M 2 be the number of mesh points of such meshes, then, for i = 1, 2, we have that where the mth element, e m := {x : r := {r ∈ R : 0 < r < R s (x l )}, which represents the spherical particle at x l , and let where N E (l) denotes the number of elements in the interval [0, R s (x l )] and ∆r on the contrary, if {x l } is a left boundary point, then and if {x l } is a right boundary point, then We define the meshes D 3h∆r as ∆r .
The families of conforming linear finite element spaces associated with these meshes are the following.
where P 0 ( e l ) is the set of polynomials of degree zero defined on e l . Let {χ l (x)} M2 l=1 be the set of nodal basis functions for V where h l denotes the length of the element e l , and the Regarding the meshes D 3h∆r , we define the finite element space V h∆r (D 3 ) as follows.
or equivalently, using the notation v .
Thus, the finite element formulation is as follows. For all t ∈ (0, T end ), the semi-discrete approximation , is solution to the following system of equations.
In this system, Note that J h also depends on t through u h , v sh and η h .

Remark 2
We must note, see (14), that v h∆r and w h∆r are elementwise constant functions in the x direction, and J h is a piecewise continuous function in x for which it makes sense to consider the approximation, h J h one readily shows, by performing the integral on D 2 , that (16) can be recast as follows: for all mesh-point Once v (l) ∆r (r, t) is known, one calculates v h∆r by the expression (14).
Notice that this equation is the finite element approximation of (11) for w ∈ H 1 r (0, R s (x l )), see [1]. Based on Remark 1 and since H 1 (D i ) ֒→ C(D i ), we introduce the following spaces which are used in the error analysis and in the application of the fixed point theorems in Section 4.
where P and Q are constants sufficiently large; for i = 1, 2, C part (D i × [0, T end ]) denotes the set of piecewise continuous functions in time and continuous in space and C * part (D 2 × [0, T end ]) denotes the set of piecewise continuous functions in both time and space. S P is the candidate pool for the concentration u and its approximate u h , whereas S Q and S * Q play the same role for the concentrations vs vmax and v sh vmax respectively. So, we make the following assumption.
A4) There exist constants P, Q and K sufficiently large such that that for almost every t ∈ [0, T end ] the following bounds hold: and for i = 1, 2,

Error analysis for the semidiscrete problem
We present in this section the error analysis for the semidiscrete potentials and concentrations. Since the development of such an analysis is long, we have split its presentation in a sequence of three subsections.
In the first one, we introduce some auxiliary results needed for the error analysis. The second subsection deals with the error estimate for the potentials. Observing that the P2D model is a nonlinear coupled system of equations, then the error for the potentials depends on the error estimates for the concentrations, the analysis of which is carried out in the last subsection.

Auxiliary results
It is well known [2] that for the finite element spaces V where it should be understood that for p = 0, ∆r [0, R] be a linear finite element space where we approximate such functions, one can prove, following the approach used to prove Lemmas 1 and 2 in [5], that when w ∈ H 2 r (0, R), We consider the interpolants I , and the elliptic projection P 1 : where λ > 0 is a constant; the error analysis for elliptic problems suggests that a good choice is λ = k 1 . By virtue of (24) it follows that there exists a constant C independent of h such that and from the well known error analysis for elliptic problems [2] Likewise, for symmetric radial functions v ∈ H q r (0, R), we define the elliptic projector P r 1 : with λ > 0; as before, a good choice now is λ = k 2 . By virtue of (25) it follows that there exists a constant C independent of ∆r such that P r 1 can be extended to functions of x and r in an L 2 -sense. Thus, for (x, r) ∈ D 2 × (0, R(·)) we define the extended projection P r 1 : We shall also consider the x-Lagrange interpolant for functions that depend on x and r, I x 0 : Since I x 0 can be viewed as an extended Lagrange interpolant I , and let f ∈ H 1 (I). There exists an arbitrarily small number ǫ and a positive constant C(ǫ) such that Proof. Since H 1 (I) ֒→ C(I), then f ∈ C(I) and so does f 2 , so for any y, z ∈ I, y < z, we have that Since exists x * ∈ I such that f 2 (x * ) = min x∈I f 2 (x), then letting y = x * it follows that Substituting this estimate the result follows. The next result is a rewording of Lemma 2.4 of [16]. Let F 1 (R) be the closure of C ∞ − functions with respect to the H 1 r (0, R)-norm and with the property that their with first derivative vanishes at r = 0.
2) There exists an arbitrarily small number ǫ and a positive (possibly large) constant C(ǫ), both depending on a, such that Proof. To prove 1) we note that for any So, any sequence {v n } that converges with respect to the H 1 r (0, R)−norm also converges with respect to the H 1 (a, R)−norm. As for the point 2), we notice that from (36) and (38) it readily follows (37).
we have the following estimates.
Proof. Recalling the expressions for J, see (2)-(3), and J h , see (20), using the assumption A4 and the bounds (22) and (23), we have that for all ( where due to the bounds (22) and (23) the constants C = C(P, Q, K). Now, by virtue of the mean value theorem there exists z ∈ (η, η h ) such that and resorting again to (22) and (23) it follows that

Hence, applying Lemma 5 yields
From this estimate it follows (40).
Furthermore, concerning the right hand side terms of (12) and (13), we introduce the operator B : V → V * , V * being the dual for V , as Hence, we can recast the equations (12) and (13) as follows. Find Φ ∈ L 2 (0, Likewise, the finite element solutions φ 1h and φ 2h that satisfy (17) and (18) respectively, can be formulated as follows. where with Remark 7 The existence and uniqueness of Φ is proven in [4] under the same kind of assumptions as A1-A4, whereas in [19] and [10] the existence is proven applying Schauder Fixed Point Theorem [7] and the uniqueness using the fact that φ 1 (x) ∈ W (D 1 ).
As for the bilinear form a and the operator B, we have the following result.
Lemma 8 Assuming that A1-A4 hold, we have that: (i) the bilinear form a is continuous, (ii) the operator B is monotone, i.e., bounded and continuous in the sense that for all Ψ ∈ V there exists a constant C such that Proof. It is easy to prove the continuity of the bilinear form a if one takes into account the regularity assumption A2. To prove (45) we note that Since for all x ∈ D 2 a 2 i 0 > 0, then by virtue of A4 we can choose a constant C(P, Q) such that for all x ∈ D 2 C(P, Q) ≤ a 2 i 0 , and by the mean value theorem sinh(βη) − sinh(β η) ≥ β(η − η), then one readily obtains To prove that B is bounded we notice that for all Φ ∈ V but |J| is bounded by virtue of A4, then using the Cauchy-Schwarz inequality it readily follows that there exists a bounded positive constant C such that so B is bounded. To prove that B is continuous, we again notice that so, arguing as in the proof of Lemma 6 we have that there exists a positive constant C such that Substituting this estimate in the above inequality and making use of the Cauchy-Schwarz inequality it follows that

Corollary 9
The discrete bilinear form a h defined in (43) is continuous. The discrete operator B h defined in (44) is monotone, bounded and continuous.
Theorem 10 For a.e. t ∈ [0, T end ], let the solution of (41), Φ(t) = (φ 1 (t), φ 2 (t)), be in H 2 (D 1 )×H 2 (D 2 ). There exists a constant C independent of h such that Proof. Setting Ψ = Ψ h in (41) and subtracting (42) yields Noting that To estimate the terms of this expression we choose Ψ h = I . For convenience, we shall split the expression for Ψ h as

Replacing this expression for Ψ h in (49) we have that
We bound the terms of (50). We start by showing that there exists a positive constant α such that the term on the left hand side satisfies To do so we note that by virtue of (47) Since φ 1 −φ 1h ∈ W (D 1 ), we can use A2 and Poincaré-Wirtinger inequality to bound a 1 (φ 1 −φ 1h , φ 1 −φ 1h ) from below as where the constant c 1 = κ 0 (1 + C P ) −1 , C P being the constant of the Poincaré-Wirtinger inequality; using again A2, we bound the term a 2 (φ 2 − φ 2h , φ 2 − φ 2h ) as Applying Young inequality we find that there exists a constant γ ∈ (0, 1) such that Now, we can choose the constants C and γ such that c 2 = c 1 + C(1 − 4 γ ) > 0, and substitute these bounds in (52) to obtain the inequality (51), where α = min(c 1 , c 2 , σ 0 , (1 − γ)C). Next, we bound the terms on the right hand side. By continuity of the bilinear form and Young inequality, we find that there exists a small positive number ǫ 1 and a constant C(ǫ 1 ) such that By virtue of (22) and (23) and the mean value theorem for the integral Applying Young inequality yields where ǫ 2 is a small positive number and C(ǫ 2 ) is a constant. Next, noting that by virtue of (22) and (23) Again, using Lemma 5 and Young inequality we obtain that there exist a small number ǫ 3 and a constant C(ǫ 3 ) such that To bound the last term on the right hand side of (50) we note that ∂φ1 ∂x L ∞ (D1) is bounded and by virtue of assumption A4, and by the mean value theorem, |(κ(u) − κ(u h ))| ≤ C |u − u h |, then it follows that Letting ǫ 1 + · · · + ǫ 4 = α/2 and noting that, see (34), the estimate (48) follows from (50)-(56).

Error estimates for the concentrations u(x, t) and v(x; r, t)
We wish to estimate u(x, t) − u h (x, t) and v(x; r, t) − v h∆r (x; r, t) in the L 2 −norm assuming that both u(x, t) and v(x; r, t) are as regular as required. Following the standard approach, we decompose u − u h as where P 1 is the elliptic projector defined in (26), and note that θ u (x, t) ∈ V (1) h (D 1 ). To carry out a decomposition of this kind for v(x; r, t) − v h∆r (x; r, t), at first we can try using the extended elliptic projector P r 1 : ∆r [0, R s (·)] defined in (32) and assume that for all t, v(x; r, t) ∈ H 1,1 r (D 2 × (0, R s (·))), then we find that for a. e. x ∈ D 2 P r 1 v(x; r, t) = ) and the function P r 1 v(x; r j , t) ∈ L 2 (D 2 ); so, in general, P r 1 v(x; r, t) is not in V h∆r (D 3 ) and, consequently, it does not make sense to use P r 1 v(x; r, t) − v h∆r (x; r, t) for such type of decomposition; however, recalling the interpolant I x 0 : ) defined in (34) and further assuming that P r 1 v(x; r j , t) ∈ H 1 (D 2 ), then it follows that P r 1 v(x; r, t) ∈ H 1,1 r (D 2 × (0, R(·))) and, therefore, we can define I x 0 P r 1 v(x; r, t) as this expression implies that I x 0 P r 1 v(x; r, t) ∈ V h∆r (D 3 ), so it makes sense to set Now, using again the extended P r 1 elliptic projector we define ρ v (x; r, t) = v(x; r, t) − P r 1 v(x; r, t) ∈ H 1,1 r (D 2 × (0, R(·))), and consequently

Then, from all these considerations we can write that
From (57) and (58) it follows that ) . The estimates for ρ u and ρ v are given in (29) and (33) respectively, i.e., then, it remains to calculate the estimates for θ u and θ v ; but before going into the details of such calculations, we present new estimates for , which depend on θ u and θ v respectively, and will be useful for the subsequent part of the analysis.

Lemma 11
Assuming that the regularity assumptions required in the estimates hold, there exist an arbitrarily small positive number ǫ and constants C and C(ǫ) independent of h and ∆r such that and Proof. To calculate the estimate (60) we set , t) (recalling the definition of θ v ) then by virtue of the definition of the L 2 -norm for functions of V 0 h (D 2 ) presented in Section 3, we have that To estimate ρ 2 v (x l ; R s (x l ), t) and θ 2 v (x l ; R s (x l ), t) we make use of Lemma 4 noting that there exists a real number a, 0 < a < R s (x l ), such that L ∞ (a,R(x l )) , thus by virtue of (37) it follows that there are a real number ǫ and positive r (0,Rs(·))) , because by approximation theory I x 0 ρ v (t) L 2 (D2,H 1 r (0,Rs(x))) ≤ C ρ v (t) L 2 (D2,H 1 r (0,Rs(x))) , this latter term being estimated according to (33). Similarly, So, putting these bounds together the result (60) follows. To calculate the estimate (61) we notice that by virtue of (40) and Theorem 10 Since u − u h = ρ u + θ u , then taking into account (29) so the result (61) follows. Next, we calculate an estimate for θ v . To this end, we obtain, based on equation (11), the integral equation for I x 0 v(x; r, t) that will be used for that purpose. Thus, for each one of the mesh points {x l } M2 l=1 of D 2h the equation (11) reads where v (l) := v(x l ; r, t) ∈ H 1 r (0, R s (x l ) and w (l) := w(x l ; r) ∈ H 1 r (0, R s (x l )). Using the nodal basis functions {χ l (x)} Rs(x) 0 v (l) (r, t)w (l) (r)r 2 dr dx, and for x ∈ e l , R s (x) = R s (x l ), then it follows that M2 l=1 D2 We proceed to formulate the equation for θ v . From (58) it follows that v h∆r = I x 0 v − (I x 0 ρ v + θ v ), then replacing this expression for v h∆r in (16) and using (32) and (64) it follows that for all w h∆r ∈ V h∆r (D 3 ), where we have made use of the following properties of the interpolant r (0,Rs(·))) + k 2 ∂θ v (t) ∂r 2 L 2 (D2,L 2 r (0,Rs(·))) ≤ λ I x 0 ρ v (t) L 2 (D2,L 2 r (0,Rs(·))) θ v (t) L 2 (D2,L 2 r (0,Rs(·))) where, θ vs (t) = θ v (x; R s (x), t) is the value of θ v on the surface of the sphere of radius R x (x) associated with the point {x} of D 2 .

Hence, by Young inequality it follows that
We bound the term R 3 . Thus, we have that Estimating the last term on the right hand side of this inequality as we did before in the proof of Lemma 11, see (62), and noting that I 0 , it readily follows that To bound R 4 we notice that J h − I , so by the triangle inequality it follows that Noting that ab L 1 (D2) ≤ ε 2 a 2 L 2 (D2) + 1 2ε b 2 L 2 (D2) and applying the same argument as we have just done to bound R 3 , we obtain that We bound the last term of this inequality as we have done for R 3 , and by virtue of assumption A3 set Letting ǫ = k 2 /8 in (61), (70) and (71), and replacing (68)-(71) in (66), the result (67) follows. Next, we proceed to calculate an estimate for θ u . Thus, subtracting (15) from (10) it readily follows that Setting w h = θ u in this equation yields . By virtue of Lemma 11 we have the following result.
We are now in a position to establish the main result of this subsection.
Considering this estimate together with those of Lemmas 12 and 13 it readily follows that r (0,Rs(·))) .

Fully discrete model
We consider now the fully discrete model based on the time stepping backward Euler scheme. This scheme has been used to discretize in time the equations of the P2D model either with finite differences [12], finite volumes [17], [11], or finite elements. For convenience, hereafter we shall use the notation a n := a(x, t n ), where n is a nonnegative integer and t n = n∆t, ∆t being a uniform time step. The formulation of the fully discrete model is as follows. Assuming that at time t n−1 , n = 1, 2, . . . , N , the solution (u n−1 where

On the existence and uniqueness of the solution of the fully discrete model
To prove that the system (75)-(78) has a unique solution, we first show that assuming (u n h , v n h∆r , v n sh ) ∈ V h (D 2 ) and the assumptions A1-A4 hold, the system (77)-(78) has a unique solution h (D 2 ); then, returning to the system (75)-(76) and applying a wellknown consequence of Brower´s fixed point theorem, which is presented as Corollary 1.1 in [8], we prove that there exists (u n h , v n h∆r ) ∈ V (1)

Lemma 15
Assuming that for all n, , and the assumptions A1-A4 hold, then the system (77)-(78) has a unique solution Proof. Looking at (77)-(78) and in order to apply Minty-Browder theorem to prove the existence of a solution, we define the functions it is worth noticing that η n h is equal to η n h when the potentials φ n 1h and φ n 1h are zero. Now, going back to Section 4.2 and using J n h , we define the operators B h : V h → V * h and A h : V h → V * h as follows: for all n = 1, 2, .., N , Notice that when Φ n h = (0, 0), B h (Φ n h ), Ψ h = 0 because J n h = J n h . Now, we can recast (77)-(78) as follows. Find Φ n h : We can prove, using the same arguments as in Lemma 8, that the operator B h is monotone, bounded and continuous satisfying an inequality as (46); since the bilinear form a h is continuous and semi-definite positive, then it follows that the operator A h is monotone, bounded and continuous satisfying an inequality as (46). In order to prove that (81) has a solution it remains to show that A h is coercive, i.e., ∀Φ n h ∈ W h (D 1 ) × V (1) h (D 2 ), there exists a positive constant α such that This can be easily done by considering the following facts: 1) A h is monotone; 2) it is easy to check, using the same arguments as in Theorem 10 to prove (51), that ∀Φ n h , Φ then taking Φ n h = (0, 0) it follows the coerciveness of A h . Hence, the Minty-Browder theorem [18] guaranties the existence of a solution Φ n h of (81). To prove the uniqueness of this solution we follow the argument put forward in [19] to prove the uniqueness of the exact solution, and assume that there two solutions Φ n h := (φ n 1h , φ n 2h ) and Φ n h := φ n 1h , φ n 2h of (77)-(78), then setting, z 1h = φ n 1h − φ n 1h and and from (78) . Setting w h = z 1h and v h = z 2h and applying the mean value theorem one readily obtains that ∂η n > 0 according to assumption A3. The first term of this expression implies that for all n, Similarly, from the second and third terms it follows that z 2h = 0, and consequently φ n 2h = φ n 2h . Hence, we have just proved that for all n there is a unique solution (φ n 1h , φ n 2h ).
Proof. We start proving the existence of u n h ∈ V h (D 1 ) as solution of (75). To this end, we write (75) h (D 1 ) is a continuous mapping defined by the relation with v sh being picked up from S * Q because we assume that v n h∆r belongs to this space; moreover, we also assume that χ h is in S P . According to Brower´s fixed point theorem, the equation On account of the assumptions v sh ∈ S * Q and χ h ∈ S P , it follows that there exists a constant Then, taking ∆t ≤ ∆t 0 < 1/C 1 , D1 F h (χ h )χ h dx is positive for χ h L 2 (D1) sufficiently large. This shows the existence of the solution u n h ∈ V h (D 1 ). Next, we prove the uniqueness. To this end, we consider that there exist X and Y ∈ V h (D 1 ).
Setting w h = X − Y and invoking the arguments of Lemmas 5 and 6 yields where the constant C 2 = C 2 (P, Q, K). Thus, taking ∆t ≤ ∆t 0 < 1/C 2 it follows that X = Y . It remains to prove the existence and uniqueness of v n h∆r , but the arguments to be used for such a proof are the same as for u n h , so we omit them.

Error estimates for the fully discrete solution
As in Section 4.3, we write for Theorem 17 Let (u n h , v n h∆r , φ n 1 , φ n 2 ) be the solution to (75)-(78). Then, under proper regularity assumptions there exists a constant C such that for ∆t small The constant C is of the form C(Γ) exp(C(k 1 , k 2 )t n , C(Γ) being another constant that depends on the exact solution (u, v, φ 1 , φ 2 ), see (87) below.
We bound the right hand side of this inequality applying the same arguments as in (66). Thus, we have that R n 1 ≤ C∆r 4 v n 2 L 2 (D2,H 2 r (0,Rs(·))) + C θ n .
Hence, we can write .
This completes the proof