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

R. Bermejo
$Dpto.Matem\'{a}ticaAplicadaalaIngenier\'{i}aIndustrialETSIIUniversidadPolit\'{e}cnicadeMadrid.$
e-mail<EMAIL_ADDRESS>
###### Abstract

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.

Keywords: P2D model, lithium-ion batteries, nonlinear, parabolic, elliptic,
finite elements, error estimates. 2010:MSC: 65M60, 35M13, 35Q99.

##  1 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 LixC6 in the negative electrode and metal oxide, such as
Li1-xMn${}_{\mathrm{2}}O_{\mathrm{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.

Figure 1: Upper panel: a cross-section of a cell along the x-direction. The
lithium ions travel from the anode to the cathode during the discharge process
and in the opposite direction during the charge process. Middel panel: the
cell model as a non-denumerable collection of solid spheres plus the
separator, there is one sphere of radius $R_{s}(x)$ at each point ${x}$ of the
electrodes. Lower panel: the domain $D_{3}$

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 $\frac{L_{x}}{L_{y}}$ and
$\frac{L_{x}}{L_{z}}=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.,
$\mu^{eff}=\mu\varepsilon^{p}$ (p=1.5), where $\mu$ is a generic transport
coefficient and $\varepsilon$ is the component volume fraction of the material
in the composite electrodes and separator.To formulate the equations of the
model we distinguish the following domains.

$\left\\{\begin{array}[]{l}D_{\mathrm{n}}=(0,L_{n}),\
D_{\mathrm{s}}=(L_{n},L_{n}+\delta)\text{,\ }D_{\mathrm{p}}=(L_{n}+\delta,L),\
L_{p}:=L-(L_{n}+\delta),\\\ \\\ D_{1}=(0,L),\ D_{2}=D_{\mathrm{n}}\cup
D_{\mathrm{p}}\text{ \ {and\ \ }}D_{3}=\cup_{x\in
D_{2}}\left\\{x\right\\}\times(0,R_{\mathrm{s}}(x)),\end{array}\right.$

where $D_{\mathrm{n}}$, $D_{\mathrm{s}}$ and $D_{\mathrm{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_{\mathrm{s}}(x)$ that represent at each $x\in
D_{2}$ the solid active particles, such that when $x\in D_{\mathrm{n}}$,
$R_{\mathrm{s}}(x)=R_{\mathrm{s}}^{-}$, and when $x\in$ $D_{\mathrm{p}}$,
$R_{\mathrm{s}}(x)=R_{\mathrm{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 $\phi_{1}(x,t)$, $x\in D_{1}$; for the
solid phase, the molar concentration of lithium $v(x;r,t)$, $x\in D_{2}$ and
$r\in\left(0;R_{\mathrm{s}}(x)\right)$, and the electric potential
$\phi_{2}(x,t)$, $x\in D_{2}$. Another important variable is the so called
molar flux of lithium ions exiting the solid particles,
$J(x,u,v,\phi_{1},\phi_{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.

##  2 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)$, $\phi_{1}(x,t)$ and $\phi_{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 $\phi_{2}$ one considers the function $H(x,t)$ given by the
expression

$H(x,t)=\left\\{\begin{array}[]{l}-\displaystyle\frac{(x-L_{n})^{2}I(t)}{2\sigma
L_{n}A},\ \ x\in D_{\mathrm{n}},\\\ \\\
\displaystyle\frac{(x-(L_{n}+\delta))^{2}I(t)}{2\sigma L_{p}A},\ \ x\in
D_{\mathrm{p}},\end{array}\right.$

where $I(t)$ denotes the applied current, $A$ is the area of the plate and
$\sigma$ is a positive coefficient defined below, and replace $\phi_{2}(x,t)$
by $\phi_{2}(x,t)+H(x,t)$; likewise, we replace the potential $\phi_{1}(x,t)$
by $\phi_{1}(x,t)+\alpha\ln u(x,t)$, with
$\alpha=\alpha(u)=\displaystyle\frac{2RT\kappa(u)}{F}(t_{+}^{0}-1)$, where
$\kappa(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 $\phi_{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.

$J=J(x,u,v_{s},T,\eta)=\left\\{\begin{array}[]{l}a_{s}i_{0}\displaystyle\left(\exp\frac{\alpha_{a}F}{RT}\left(\eta-\frac{R_{SEI}}{a_{s}}J\right)\right.\\\
\\\
-\displaystyle\left.\exp\frac{-\alpha_{c}F}{RT}\left(\eta-\frac{R_{SEI}}{a_{s}}J\right)\right)\
\ \text{if \ }x\in D_{n}\cup D_{p},\\\ \\\ 0\ \ \text{if \ }x\in
D_{s}.\end{array}\right.$ (1)

In this expression, $v_{s}=v(x;R_{s}(x),t)$ denotes the lithium concentration
on the surface of the active particles;
$a_{s}=a_{s}(x)=\displaystyle\frac{3\varepsilon_{s}(x)}{R_{s}(x)}$ is the
active area per electrode unit volume; $\varepsilon_{s}(x)$ denotes the volume
fraction of the active material, $\varepsilon_{s}(x)=\varepsilon_{s}^{-}>0$
for $x\in D_{n}$ and $\varepsilon_{s}(x)=\varepsilon_{s}^{+}>0$ for $x\in
D_{p}$; $\alpha_{a}\in(0,1)$ and $\alpha_{c}\in(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$.

$\eta=\eta(x,\phi_{1},\phi_{2},U)=\left\\{\begin{array}[]{l}\phi_{2}(x,t)-\phi_{1}(x,t)-U(x,v_{s})\text{
if\ }x\in D_{n}\cup D_{p},\\\ \\\ 0\text{ {if} }x\in D_{s},\end{array}\right.$

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.,

$i_{0}=i_{0}(u,v_{s})=ku^{\alpha_{a}}(v_{\max}-v_{s})^{\alpha_{a}}v_{s}^{\alpha_{c}}\text{
if }x\in D_{n}\cup L_{n},$ (2)

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

$v_{\max}=v_{\mathrm{\max}}(x)=\left\\{\begin{array}[]{l}v_{\mathrm{\max}}^{-}\text{
{if} }x\in D_{n},\\\ v_{\mathrm{\max}}^{+}\text{ {if} }x\in
D_{p},\end{array}\right.$

the coefficient $k$ represents the kinetic rate constant,

$k=k(x)=\left\\{\begin{array}[]{l}k^{-}\text{ {if} }x\in D_{n},\\\ k^{+}\text{
{if} }x\in D_{p}.\end{array}\right.$

Considering the above mentioned changes of variable and taking the transfer
coefficients $\alpha_{a}$ and $\alpha_{c}$ equal to $0.5$, as many engineering
papers do, the expression for the reaction current that we use in the paper is

$J=J(x,u,v_{s},\phi_{1},\phi_{2},\overline{U})=\left\\{\begin{array}[]{l}a_{2}(x)i_{0}\sinh\left(\beta\eta\right),\
\forall x\in D_{2},\\\ \\\ 0\ \mathrm{for\ }x\notin D_{2},\end{array}\right.$
(3)

where $\beta=\displaystyle\frac{F}{2RT}$;
$a_{2}(x)=3\varepsilon_{s}(x)/R_{s}(x)$; and

$\eta=\phi_{2}-\phi_{1}-\alpha\ln u-\overline{U},$ (4)

$\overline{U}=\overline{U}(x,t,v_{s})=U(v_{s})-H(x,t)$. Noting that the
boundaries $\partial D_{1}$ and $\partial D_{2}$ of the domains $D_{1}$ and
$D_{2}$ are $\partial D_{1}:=\left\\{0,L\right\\}$ and $\partial
D_{2}:=\left\\{0,L_{n},L_{n}+\delta,L\right\\}$, we formulate the equations of
the model as follows.

Concentration $u(x,t)$ in the electrolyte phase.

$\left\\{\begin{array}[]{l}\displaystyle\frac{\partial u}{\partial
t}-\displaystyle\frac{\partial}{\partial x}(k_{1}\frac{\partial u}{\partial
x})=a_{1}(x)J\text{ \ {in} \ }D_{1}\times(0,T_{\text{{end}}}),\\\ \\\
\displaystyle\frac{\partial u}{\partial x}\mid_{\partial
D_{1}\times(0,T_{\text{{end}}})}=0,\ u(x,0)=u^{0}(x)\ \text{{in}
}D_{1}.\end{array}\right.$ (5)

Concentration $v(x;\,r,t)$ in the solid phase. For almost every $x\in D_{2}$,

$\left\\{\begin{array}[]{l}\displaystyle\frac{\partial v}{\partial
t}-\frac{k_{2}}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial
v}{\partial r}\right)v=0\text{\ \ {in}\ \ }D_{3}\times(0,T_{\text{{end}}}),\\\
\\\ \displaystyle\frac{\partial v}{\partial r}\mid_{r=0}=0,\
k_{\mathrm{2}}\displaystyle\frac{\partial v}{\partial
r}\mid_{r=R_{\mathrm{s}}(x)}=\displaystyle\frac{-J}{a_{2}(x)F},\
v(x;r,0)=v^{0}(x;r)\ \ \text{{in} \ }D_{3}.\end{array}\right.$ (6)

Electrolyte potential $\phi_{\mathrm{1}}(x,t)$.

$\left\\{\begin{array}[]{l}-\displaystyle\frac{\partial}{\partial
x}(\kappa(u)\frac{\partial\phi_{\mathrm{1}}}{\partial x})=J\text{ \ {in} \
}D_{1}\times(0,T_{\text{end}}),\\\ \\\
\displaystyle\frac{\partial\phi_{\mathrm{1}}}{\partial x}\mid_{\partial
D_{1}\times(0,T_{\text{end}})}=0,\\\ \\\
\int_{D_{1}}\phi_{\mathrm{1}}(x,t)dx=0.\end{array}\right.$ (7)

Solid phase potential $\phi_{\mathrm{2}}(x,t)$.

$\left\\{\begin{array}[]{l}\displaystyle\frac{\partial}{\partial
x}\left(\sigma\frac{\partial\phi_{\mathrm{2}}}{\partial x}\right)=J+g\text{ \
{in} \ }D_{2}\times(0,T_{\text{{end}}}),\\\ \\\
\sigma\displaystyle\frac{\partial\phi_{\mathrm{2}}}{\partial x}\mid_{\partial
D_{2}\times(0,T_{\text{{end}}})}=0,\\\ \\\
g(x,t)=\left\\{\begin{array}[]{r}\displaystyle\frac{-I(t)}{L_{n}A},\ x\in
D_{\textrm{n}},\\\ \\\ \displaystyle\frac{I(t)}{L_{p}A},\ x\in
D_{\textrm{p}},\end{array}\right.\end{array}\right.$ (8)

where $a_{1}(x)=\displaystyle\frac{1-t_{+}^{0}}{3\varepsilon_{s}(x)F}$. In
these equations, $k_{1}(x)>0$ and $k_{2}(x)>0$ represent effective diffusion
coefficients in the electrolyte and solid phases respectively, and $\sigma(x)$
denotes the effective electric conductivity in the solid phase. The functions
$a_{1}(x),\ a_{2}(x)\ $, $\sigma(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\in(0,T_{\text{{end}}})$,
$J(x,u,v_{\mathrm{s}},\phi_{\mathrm{1}},\phi_{\mathrm{2}})$ satisfies the
algebraic conditions

$\left\\{\begin{array}[]{l}\displaystyle\int_{D_{1}}Jdx=\int_{D_{2}}Jdx=0,\\\
\\\ \displaystyle\int_{D_{\textrm{n}}}Jdx=I(t),\
\int_{D_{\textrm{p}}}Jdx=-I(t).\end{array}\right.$ (9)

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\in\left[0,T_{\mathrm{end}}\right]$

$\int_{D_{1}}u(x,t)dx=\int_{D_{1}}u(x,0)dx,$

and

$\int_{D_{2}}\int_{0}^{R_{s}(x)}v(x;r,t)r^{2}drdx=\int_{D_{2}}\int_{0}^{R_{s}(x)}v(x;r,0)r^{2}drdx.$

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\geq 0$ and $x\in D_{1}$, $u(x,t)>0$, similarly, for $(x,r)\in D_{3}$,
$0<v(x;r,t)<v_{\max}$.

Let $D$ denote a generic open bounded domain in $\mathbb{R}$; hereafter, the
closure of a domain $D$ is denoted $\overline{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\leq p\leq\infty$; the spaces of measurable radial
functions [4]

$H_{r}^{q}(0,R):=\left\\{v:(0,R)\rightarrow\mathbb{R}\text{:{\
}}\left\|v\right\|_{H_{r}^{q}(0,R)}^{2}=\sum_{j=0}^{q}\int_{0}^{R}\left(\frac{d^{j}v}{dr^{j}}\right)^{2}r^{2}dr<\infty\right\\},$

$q$ being a nonnegative integer, when $q=0$ we set
$H_{r}^{0}(0,R):=L_{r}^{2}(0,R)$; also, for $p$ being a nonnegative integer,
the normed spaces of measurable functions

$H^{p}(D_{2};H_{r}^{q}(0,R_{s}(\cdot))):=\left\\{v:D_{2}\rightarrow
H_{r}^{q}(0,R_{s}(\cdot)):\left\|v\right\|_{H^{p}(D_{2};H_{r}^{q}(0,R_{s}(\cdot)))}<\infty\right\\},$

where $\left\|v\right\|_{H^{p}(D_{2};H_{r}^{q}(0,R_{s}(\cdot)))}^{2}=$
$\displaystyle\sum_{j=0}^{p}\int_{D_{2}}\left\|\frac{\partial^{j}v(x;\cdot)}{\partial
x^{j}}\right\|_{H_{r}^{q}(0,R_{s}(x))}^{2}dx$; and the spaces

$H_{\ \
r}^{p,q}(D_{2}\times(0,R_{s}(\cdot)))=H^{p}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))\cap
L^{2}(D_{2},H_{r}^{q}(0,R_{s}(\cdot)))$

with norm

$\left\|u\right\|_{H_{\ \
r}^{p,q}(D_{2}\times(0,R_{s}(\cdot)))}^{2}:=\left\|u\right\|_{H^{p}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\left\|u\right\|_{L^{2}(D_{2},H_{r}^{q}(0,R_{s}(\cdot)))}^{2}.$

Notice that $H_{\ \
r}^{0,0}(D_{2}\times(0,R_{s}(\cdot)))=L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot))).$

Since the variables of the model depend on time, then we also introduce the
normed spaces $L^{p}(0,t;X)$, where $1\leq p\leq\infty$, and $(X,\
\left\|\cdot\right\|_{X})$ being a real Banach space.

$L^{p}(0,t;X):=\left\\{v:(0,t)\rightarrow X\text{ {strongly measurable such
that\ \ }}\left\|v\right\|_{L^{p}(0,t;X)}<\infty\right\\},$

with
$\left\|v\right\|_{L^{p}(0,t;X)}=\displaystyle\left(\int_{0}^{t}\left\|v(\tau)\right\|_{X}^{p}d\tau\right)^{1/p}$
when $1\leq p<\infty$, and for $p=\infty$,
$\left\|v\right\|_{L^{\infty}(0,t;X)}=ess\
\sup_{0<\tau<t}\left\|v(\tau)\right\|_{X}$. Other spaces used in the paper are
$W(D_{1}):=\left\\{v\in H^{1}(D_{1}):\displaystyle\int_{D_{1}}vdx=0\right\\}$,
which is a closed subspace of $H^{1}(D_{1})$ where the potential
$\phi_{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)$.

Next, we introduce the following regularity assumptions on the data and the
molar flux $J$ [4].

A1)

$u^{0}\in H^{1}(D),\ u^{0}>0,\ v^{0}\in C(\overline{D}_{3}),\
0<v^{0}<v_{\max},\ I(t)\in C_{part}([0,T^{\ast}]),0<T_{end}\leq
T^{\ast}<\infty,$

where $C_{part}$ denotes the set of piecewise continuous functions, i.e.,

$C_{\mathrm{part}}([a,b])=\left\\{g:[a,b]\rightarrow\mathbb{R}:\exists
a=t_{0}<t_{1}<\cdots t_{N}=b\text{ such that }g\in
C\left([t_{i-1},t_{i}]\right)\right\\}.$

A2) For $0<a<c<+\infty$, $k_{0}$ and $\sigma_{0}$ positive constants,

$k_{1}\in L^{\infty}(D_{1}),\ k_{1}\geq k_{0}>0,\ k_{2}\in[a,c],\ \kappa\in
C^{2}\left((0,+\infty)\right),\ \sigma\in L^{\infty}(D_{2}),\
\kappa\geq\kappa_{0}>0,\ \sigma\geq\sigma_{0}>0.$

Moreover, $\underline{k_{1}}:=\inf_{x\in{D_{1}}}{k_{1}(x)}$ and
$\underline{k_{2}}:=\min_{x\in{D_{2}}}{k_{2}(x)}$.

A3) For all $(x,u,v_{s},\eta)\in
D_{2}\times(0,+\infty)\times(0,v_{s,\mathrm{\max}})\times\mathbb{R}$,

$J\in
C^{2}(D_{2}\times(0,+\infty)\times(0,v_{s,\mathrm{\max}})\times\mathbb{R}),\
\displaystyle\frac{\partial J}{\partial\eta}>0.$

A weak formulation to (5)-(8) is the following. Find

$\left\\{\begin{array}[]{l}u\in{L}^{2}\left(0,T_{\mathrm{end}};\
H^{1}(D_{1}\right)),\ \displaystyle\frac{du}{dt}\in L^{2}(0,T_{\mathrm{end}};\
H^{1}(D_{1}))^{\ast},\\\ \\\ v{\in
L}^{2}\left(0,T_{\mathrm{end}};L^{2}(D_{2,}\
H_{r}^{1}(0,R_{s}(\cdot))\right)),\
\displaystyle\frac{dv}{dt}\in{L}^{2}\left(0,T_{\mathrm{end}};L^{2}(D_{2,}\
H_{r}^{1}(0,R_{s}(\cdot))\right))^{\ast},\\\ \\\ \phi_{\mathrm{1}}\in
L^{2}(0,T_{\mathrm{end}};\ W(D_{1}))\text{ {and} }\phi_{\mathrm{2}}\in
L^{2}(0,T_{\mathrm{end}};\ H^{1}(D_{2})),\end{array}\right.$

such that

$\int_{D_{1}}\frac{\partial{u}}{\partial{t}}wdx+\int_{D_{1}}k_{1}\frac{\partial
u}{\partial x}\frac{dw}{dx}dx=\int_{D_{1}}a_{1}Jwdx\ \forall w\in
H^{1}(D_{1});$ (10)

for a.e.$x\in D_{2}$ and for all $w\in H_{r}^{1}(0,R_{s}(x))$ radially
symmetric

$\int_{0}^{R_{s}(x)}\frac{\partial{v}}{\partial{t}}wr^{2}dr+\int_{0}^{R_{s}(x)}k_{2}\frac{\partial
v}{\partial r}\frac{\partial w}{\partial
r}r^{2}dr=\frac{-R_{\mathrm{s}}^{2}(x)Jw(R_{s}(x))}{a_{2}(x)F};$ (11)

$\int_{D_{1}}\kappa(u)\frac{\partial\phi_{1}}{\partial
x}\frac{dw}{dx}dx=\int_{D_{1}}Jwdx\text{ \ }\forall w\in H^{1}(D_{1});$ (12)

and

$\int_{D_{2}}\sigma\frac{\partial\phi_{2}}{\partial
x}\frac{dw}{dx}dx=-\int_{D_{2}}\left(J+g\right)wdx,\forall w\in H^{1}(D_{2}),$
(13)

where $L^{2}(0,T_{\mathrm{end}};\ H^{1}(D_{1}))^{\ast}$ and
${L}^{2}\left(0,T_{\mathrm{end}};L^{2}(D_{2,}\
H_{r}^{1}(0,R_{s}(\cdot))\right))^{\ast}$ denote the respective dual spaces of
$L^{2}(0,T_{\mathrm{end}};\ H^{1}(D_{1}))$ and
${L}^{2}\left(0,T_{\mathrm{end}};L^{2}(D_{2,}\
H_{r}^{1}(0,R_{s}(\cdot))\right))$.

###### 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}<\cdots<t_{N}$ of
$[0,T_{\mathrm{end}}]$, $T_{\mathrm{end}}$ being small enough, that there is a
unique weak solution $\left(u,v,\phi_{1},\phi_{2}\right)$ in each interval
$[t_{n},t_{n+1}]$, such that $\left(u,v,\phi_{1},\phi_{2}\right)\in
C([t_{n},t_{n+1}];K_{Z})$, where $K_{Z}:=H^{1}(D_{1})\times L^{2}(D_{2,}\
H_{r}^{1}(0,R_{s}(\cdot)))\times W(D_{1})\times H^{1}(D_{2})$.
$\left(u(t_{n}),v(t_{n})\right)$ being the initial condition in such an
interval. Also, Kr$\ddot{\text{o}}$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)$, $\phi_{\mathrm{1}}(x,t)$
and $\phi_{\mathrm{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}+\delta$ 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}\subset D_{1}$, we choose the family of meshes $D_{2h}$ as a subset of
$D_{1h}$. Let $NE_{1}$ and $NE_{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

$D_{ih}=\left\\{e_{m}\right\\}_{m=1}^{NE_{i}}\ \text{{and}}\
\overline{D}_{i}=\cup_{m=1}^{NE_{i}}e_{m},$

where the $m$th element, $e_{m}:=\left\\{x:x_{1}^{m}\leq x\leq
x_{2}^{m}\right\\}$, and $h_{m}:=x_{2}^{m}-x_{1}^{m}$ is the length of the
element $e_{m}$; the points $\ x_{1}^{m}$ $x_{2}^{m}$ are denoted element
nodes. We set $h=\max_{m}h_{m}$, and $\gamma=h^{-1}\min_{m}h_{m}$. The
parameter $\gamma$ is a measure of the uniformity of the meshes. The
collection of all the element nodes defines the set of nodes,
$\left\\{x_{l}\right\\}_{l=1}^{M_{i}}$, of the mesh $D_{ih}$. To construct the
family of meshes $\widehat{D}_{3h\Delta r}$ on $D_{3}$, we recall that
$D_{3}=\cup_{x\in D_{2}}\left\\{x\right\\}\times(0,R_{s}(x))$, where
$R_{s}(x)$ is the radius of the solid spherical particle associated with the
point $\\{x\\}$. Thus, for each mesh point $\left\\{x_{l}\right\\}\in D_{2h}$
we define the radial vertical domain
$D_{r}^{(l)}:=\left\\{r\in\mathbb{R}:0<r<R_{s}(x_{l})\right\\}$, which
represents the spherical particle at $x_{l}$, and let

$D_{\Delta r}^{(l)}=\left\\{e_{k}^{(l)}\right\\}_{k=1}^{NE^{(l)}}\text{ {such
that }}\overline{D}_{r}^{(l)}=\cup_{k=1}^{NE^{(l)}}e_{k}^{(l)},$

where $NE^{(l)}$ denotes the number of elements in the interval
$[0,R_{s}(x_{l})]$ and $\Delta r_{k}^{(l)}$ is the width of the element
$e_{k}^{(l)}:=\\{r:r_{1}^{(l)k}\leq r\leq r_{2}^{(l)k}\\}$, we set $\Delta
r=\max_{l}(\max_{k}\Delta r_{k}^{(l)})$ and $\gamma_{r}=\Delta
r^{-1}\min_{l}(\min_{k}\Delta r_{k}^{(l)}).$ The set of mesh points in each
mesh $D_{\Delta r}^{(l)}$ is denoted
$\left\\{r_{j}^{(l)}\right\\}_{j=1}^{M^{(l)}}$. Furthermore, let
$\left\\{\widehat{e}_{l}\right\\}_{l=1}^{M_{2}}$ be the collection of
nonconforming elements of the mesh $\widehat{D}_{2h}$ which are associated
with the nodes $\left\\{x_{l}\right\\}$, they are defined as follows: if
$\left\\{x_{l}\right\\}$ is not a boundary point, then

$\widehat{e}_{l}:=\left\\{x\in D_{2}:x_{l}-\frac{h_{l-1}}{2}\leq
x<x_{l}+\frac{h_{l}}{2}\right\\};$

on the contrary, if $\left\\{x_{l}\right\\}$ is a left boundary point, then

$\widehat{e}_{l}:=\left\\{x\in\overline{D}_{2}:x_{l}\leq
x<x_{l}+\frac{h_{l}}{2}\right\\},$

and if $\left\\{x_{l}\right\\}$ is a right boundary point, then

$\widehat{e}_{l}:=\left\\{x\in\overline{D}_{2}:x_{l}-\frac{h_{l-1}}{2}\leq
x\leq x_{l}\right\\}.$

We define the meshes $\widehat{D}_{3h\Delta r}$ as

$\widehat{D}_{3h\Delta r}:=\left\\{\widehat{e}_{l}\times\overline{D}_{\Delta
r}^{(l)}\right\\}_{l=1}^{M_{2}}\ \text{{such that
}}\overline{D}_{3}:=\cup_{l=1}^{M_{2}}\widehat{e}_{l}\times\overline{D}_{\Delta
r}^{(l)}.$

Figure 2: Panel (a): the mesh for the domain $D_{1}$, which includes the
negative electrode (anode in the figure) $D_{a}$, the separator $D_{s}$ and
the positive electrode (cathode in the figure) $D_{c}$. Panel (b): the mesh
for the domain $D_{2}=D_{a}\cup D_{c}$. Panel (c): the mesh of nonconforming
elements for the domain $D_{2}$. Panel (d): the mesh for the domain $D_{3}$.

The families of conforming linear finite element spaces associated with these
meshes are the following. For $i=1,2,$

$V_{h}^{(1)}(\overline{D}_{i}):=\left\\{v_{h}\in C(\overline{D}_{i}):\forall
e_{m}\in D_{ih},\ v_{h}\mid_{e_{m}}\in P_{1}(e_{m})\right\\},$

where $P_{1}(e_{m})$ denotes the set of linear polynomials defined on $e_{m}$.
Let $\left\\{\psi_{l}(x)\right\\}_{l=1}^{M_{i}}$ be the set of nodal basis
functions for the space $V_{h}^{(1)}(\overline{D}_{i})$, then any function
$v_{h}\in V_{h}^{(1)}(\overline{D}_{i})$ can be written as

$v_{h}(x)=\sum_{l=1}^{M_{i}}V_{l}\psi_{l}(x),\ \mathrm{where\
}V_{l}=v_{h}(x_{l})\text{.}$

Note that $V_{h}^{(1)}(\overline{D}_{i})\subset H^{1}(D_{i})$. The
nonconforming finite element space associated with the mesh $\widehat{D}_{2h}$
is defined as

$V_{h}^{(0)}(\overline{D}_{2}):=\left\\{v_{h}\in
L^{2}(D_{2}):\forall\widehat{e}_{l}\in\widehat{D}_{2h},\
v_{h}\mid_{\widehat{e}_{l}}\in P_{0}(\widehat{e}_{l})\right\\},$

where $P_{0}(\widehat{e}_{l})$ is the set of polynomials of degree zero
defined on $\widehat{e}_{l}$. Let $\left\\{\chi_{l}(x)\right\\}_{l=1}^{M_{2}}$
be the set of nodal basis functions for $V_{h}^{(0)}(\overline{D}_{2})$,

$\chi_{l}(x)=\left\\{\begin{array}[]{r}1\text{ \ {if} \
}x\in\widehat{e}_{l},\\\ \\\ 0\ \ \text{{otherwise}},\end{array}\right.$

then any function $v_{h}(x)\in V_{h}^{(0)}(\overline{D}_{2})$ is expressed as

$v_{h}(x)=\sum_{l=1}^{M_{2}}V_{l}\chi_{l}(x),\ \mathrm{where\
}V_{l}=v_{h}(x_{l}).$

It is worth remarking that for $1\leq p<\infty$, the $L^{p}$-norm of
$v_{h}(x)\in V_{h}^{(0)}(\overline{D}_{2})$ is given as
$\left\|v_{h}\right\|_{L^{p}(D_{2})}^{p}=\sum_{l=1}^{M_{2}}\widehat{h}_{l}V_{l}^{p}$,
where $\widehat{h}_{l}$ denotes the length of the element $\widehat{e}_{l}$,
and the $L^{2}$-inner product of $v_{h},w_{h}\in
V_{h}^{(0)}(\overline{D}_{2})$,
$\displaystyle\int_{D_{2}}v_{h}w_{h}dx=\sum_{l=1}^{M_{2}}\widehat{h}_{l}V_{l}W_{l}$.
Next, we introduce the finite element space $V_{\Delta
r}^{(1)}(\overline{D}_{r}^{(l)})$. For $1\leq l\leq M_{2}$,

$V_{\Delta r}^{(1)}(\overline{D}_{r}^{(l)}):=\left\\{v_{\Delta r}\in
C(\overline{D}_{r}^{(l)}):\forall e_{k}^{(l)}\in D_{\Delta r}^{(l)},\
v_{\Delta r}^{(l)}(r)\mid_{e_{k}^{(l)}}\in P_{1}(e_{k}^{(l)})\right\\}.$

So, if $\left\\{\alpha_{j}^{(l)}(r)\right\\}_{j=1}^{M^{(l)}}$ denotes the set
of nodal basis of $V_{\Delta r}(\overline{D}_{r}^{(l)})$, any function
$v_{\Delta r}^{(l)}\in V_{\Delta r}^{(1)}(\overline{D}_{r}^{(l)})\subset
H_{r}^{1}(0,R_{s}(x_{l}))$ can be written as

$v_{\Delta r}^{(l)}(r)=\sum_{j=1}^{M^{(l)}}V_{j}^{(l)}\alpha_{j}^{(l)}(r),\
\mathrm{where\ }V_{j}^{(l)}=v_{h}^{(l)}(r_{j})\text{.}$

Regarding the meshes $\widehat{D}_{3h\Delta r}$, we define the finite element
space $V_{h\Delta r}(\overline{D}_{3})$ as follows. For $1\leq l\leq M_{2}$
and $1\leq k\leq NE^{(l)}$

$V_{h\Delta r}(\overline{D}_{3}):=\displaystyle\left\\{{v}_{h\Delta r}{\in
H}_{\ \ r}^{0,1}{(D_{2}\times(0,R_{s})(\cdot)):v_{h\Delta
r}(x;r)\mid_{\widehat{e}_{l}\times e_{k}^{(l)}}\in
P_{0}(\widehat{e}_{l})\otimes P_{1}(e_{k}^{(l)})}\right\\},$

noting that when $x\in\widehat{e}_{l}$, $R_{s}(x)=R_{s}(x_{l})$. Hence, any
function $v_{h\Delta r}(x;r)\in V_{h\Delta r}(\overline{D}_{3})$ is of the
form

$v_{h\Delta
r}(x;r)=\sum_{l=1}^{M_{2}}\sum_{j=1}^{M^{(l)}}V_{lj}\chi_{l}(x)\alpha_{j}^{(l)}(r),\
\mathrm{where\ }V_{lj}=v_{h\Delta r}(~{}x_{l},r_{j}),$

or equivalently, using the notation $v_{\Delta r}^{(l)}(r)$ to denote
$v_{h\Delta r}(x_{l};r)$, we can write

$v_{h\Delta r}(x;r)=\sum_{l=1}^{M_{2}}v_{\Delta r}^{(l)}(r)\chi_{l}(x).$ (14)

. The function $v_{sh}(x):=v_{h\Delta r}(x;R_{s}(x))$ is given by the
expression

$v_{sh}(x;R_{s}(x))=\sum_{l=1}^{M_{2}}V_{lM^{(l)}}\chi_{l}(x),$

so that $v_{sh}(x)\in V_{h}^{(0)}(\overline{D}_{2})$. We calculate
$\phi_{1h}(x,t)$, which is the approximation to $\phi_{1}(x,t)$, in the finite
dimensional space

$W_{h}(\overline{D}_{1}):=\left\\{v_{h}\in
V_{h}^{(1)}(\overline{D}_{1}):\int_{D_{1}}v_{h}dx=0\right\\},\
W_{h}(\overline{D}_{1})\subset W(D_{1}).$

Thus, the finite element formulation is as follows. For all
$t\in(0,T_{\mathrm{end}})$, the semi-discrete approximation
$(u_{h}(t),v_{h\Delta r}(t),\phi_{1h}(t),\phi_{2h}(t))\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})\times
W_{h}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2}),$
$(u_{h}^{0},v_{h}^{0})\in V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta
r}(\overline{D}_{3})$, is solution to the following system of equations.

$\int_{D_{1}}\frac{\partial
u{{}_{h}}}{\partial{t}}w_{h}dx+\int_{D_{1}}k_{1}\frac{\partial u_{h}}{\partial
x}\frac{dw_{h}}{dx}dx=\int_{D_{1}}a_{1}J_{h}w_{h}dx\ \forall w_{h}\in
V_{h}^{(1)}(\overline{D}_{1}).$ (15)

$\left\\{\begin{array}[]{l}\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\displaystyle\frac{\partial
v{{}_{h\Delta r}}}{\partial{t}}w_{h\Delta
r}r^{2}drdx+\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}k_{2}\displaystyle\frac{\partial
v{{}_{h\Delta r}}}{\partial{r}}\frac{\partial w_{h\Delta r}}{\partial
r}r^{2}drdx\\\ \\\
=-\displaystyle\int_{D_{2}}\frac{R_{s}^{2}(x)J_{h}w_{h\Delta
r}(x,R_{s}(x))}{a_{2}(x)F}dx\ \ \forall w_{h\Delta r}\in V_{h\Delta
r}(\overline{D}_{3}).\end{array}\right.$ (16)

$\int_{D_{1}}\kappa(u_{h})\frac{\partial\phi_{1h}}{\partial
x}\frac{dw_{h}}{dx}dx=\int_{D_{1}}J_{h}w_{h}dx\ \ \forall w_{h}\in
V_{h}^{(1)}(\overline{D}_{1}).$ (17)

$\int_{D_{2}}\sigma\frac{\partial\phi_{2h}}{\partial
x}\frac{dw_{h}}{dx}dx=-\int_{D_{2}}\left(J_{h}+g\right)w_{h}dx\ \ \forall
w_{h}\in V_{h}^{(1)}(\overline{D}_{2}).$ (18)

$\int_{D_{1}}J_{h}dx=\int_{D_{2}}J_{h}dx=0,\text{ with
}\int_{D_{\mathrm{n}}}J_{h}dx=I(t)=-\int_{D_{\textrm{p}}}J_{h}dx.$ (19)

In this system,

$\begin{array}[]{l}J_{h}:=J(x,u_{{}_{h}},v_{sh},\eta_{h})=a_{2}(x)i_{0h}\sinh\left(\beta\eta_{h}\right),\
i_{0h}:=i_{0}(u_{h},v_{sh})\\\ \\\ \eta_{h}:=\phi_{1h}-\phi_{2h}-\alpha_{h}\ln
u_{h}-\overline{U}_{h}(v_{sh}),\ \alpha_{h}=\alpha(u_{h}).\end{array}$ (20)

Note that $J_{h}$ also depends on $t$ through $u_{h}$, $v_{sh}$ and
$\eta_{h}$.

###### Remark 2

We must note, see (14), that $v_{h\Delta r}$ and $w_{h\Delta 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, $I_{h}^{(0)}J_{h}\in V_{h}^{(0)}(\overline{D}_{2})$, see in
Subsection 4.1 the definition of the interpolant $I_{h}^{(0)}$. Then,
approximating $J_{h}$ by $I_{h}^{(0)}J_{h}$ one readily shows, by performing
the integral on $D_{2}$, that (16) can be recast as follows: for all mesh-
point $\left\\{x_{l}\right\\}\in D_{2h}$, calculate $v_{\Delta r}^{(l)}(r,t)$
$\in V_{\Delta r}^{(1)}(\overline{D}_{r}^{(l)})$ such that

$\begin{array}[]{r}\displaystyle\int_{0}^{R_{s}(x_{l})}\frac{\partial
v_{\Delta r}^{(l)}(r,t)}{\partial t}w_{\Delta
r}^{(l)}(r)r^{2}dr+\int_{0}^{R_{s}(x_{l})}\displaystyle k_{2}\frac{\partial
v_{\Delta r}^{(l)}(r,t)}{\partial r}\frac{\partial w_{\Delta
r}^{(l)}(r)}{\partial r}r^{2}dr\\\ \\\
=-R_{s}^{2}(x_{l})(a_{2}(x_{l})F)^{-1}J_{h}\mid_{x_{l}}w_{\Delta
r}^{(l)}(R_{s}(x_{l}))\ \ \forall w_{\Delta r}^{(l)}\in V_{\Delta
r}^{(1)}(\overline{D}_{r}^{(l)}).\end{array}$ (21)

Once $v_{\Delta r}^{(l)}(r,t)$ is known, one calculates $v_{h\Delta r}$ by the
expression (14).

Notice that this equation is the finite element approximation of (11) for
$w\in H_{r}^{1}(0,R_{s}(x_{l}))$, see [1].

Based on Remark 1 and since $H^{1}(D_{i})\hookrightarrow C(\overline{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.

$\begin{array}[]{l}S_{P}:=\left\\{u\in
C_{\mathrm{part}}(\overline{D}_{1}\times[0,T_{\mathrm{end}}]):\displaystyle\frac{1}{P}\leq
u\leq P\right\\},\\\ \\\ S_{Q}:=\left\\{w\in
C_{\mathrm{part}}(\overline{D}_{2}\times[0,T_{\mathrm{end}}]):\displaystyle\frac{1}{Q}\leq
w\leq 1-\frac{1}{Q}\right\\},\\\ \\\ S_{Q}^{\ast}:=\left\\{w\in
C_{\mathrm{part}}^{\ast}(\overline{D}_{2}\times[0,T_{\mathrm{end}}]):\frac{1}{Q}\leq
w\leq 1-\frac{1}{Q}\right\\},\end{array}$

where $P$ and $Q$ are constants sufficiently large; for $i=1,2$,
$C_{\mathrm{part}}(\overline{D}_{i}\times[0,T_{\mathrm{end}}])$ denotes the
set of piecewise continuous functions in time and continuous in space and
$C_{\mathrm{part}}^{\ast}(\overline{D}_{2}\times[0,T_{\mathrm{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}^{\ast}$ play the same role for the concentrations
$\frac{v_{s}}{v_{\max}}$ and $\frac{v_{sh}}{v_{\max}}$ respectively. So, we
make the following assumption.

A4) There exist constants $P,\ Q$ and $K$ sufficiently large such that that
for almost every $t\in[0,T_{\mathrm{end}}]$ the following bounds hold:

$\frac{1}{P}\leq u(t),\ u_{h}(t)\leq P,\text{ \ }\
\frac{1}{Q}\leq\frac{v_{s}(t)}{v_{\max}},\
\frac{v_{sh}(t)}{v_{\max}}\leq\left(1-\frac{1}{Q}\right),$ (22)

and for $i=1,2$,

$\left\|\frac{\partial\phi_{1}(t)}{\partial x}\right\|_{L^{\infty}(D_{1})},\
\left\|\phi_{i}(t)\right\|_{H^{1}(D_{i})},\
\left\|\phi_{i}(t)\right\|_{L^{\infty}(D_{i})},\
\left\|\phi_{ih}(t)\right\|_{L^{\infty}(D_{i})}\leq K.$ (23)

##  4 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.

###  4.1 Auxiliary results

It is well known [2] that for the finite element spaces
$V_{h}^{(p)}(\overline{D}_{i})$ ($i=1,2;p=0,1$) the following approximation
property holds. For $1\leq s\leq p+1$,

$\inf_{v_{h}\in
V_{h}^{(p)}(D_{i})}\left\\{\left\|v-v_{h}\right\|_{L^{2}(D_{i})}+h\left\|\frac{d^{p}\left(v-v_{h}\right)}{dx^{p}}\right\|_{L^{2}(D_{i})}\right\\}\leq
Ch^{s}\left|v\right|_{H^{s}(D_{i})},$ (24)

where it should be understood that for $p=0$,
$\displaystyle\frac{d^{p}\left(v-v_{h}\right)}{dx^{p}}=v-v_{h}$. For symmetric
radial functions defined in the interval $[0,R]$, let $V_{\Delta
r}^{(1)}[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\in H_{r}^{2}(0,R)$,

$\inf_{w_{\Delta r}\in V_{\Delta r}^{(1)}[0,R]}\left\\{\left\|w-w_{\Delta
r}\right\|_{L_{r}^{2}(0,R)}+\Delta r\left\|\frac{d\left(w-w_{\Delta
r}\right)}{dr}\right\|_{L_{r}^{2}(0,R)}\right\\}\leq C\Delta
r^{2}\left|w\right|_{H_{r}^{2}(0,R)}.$ (25)

We consider the interpolants $I_{h}^{(1)}:C(\overline{D}_{i})\rightarrow$
$V_{h}^{(1)}(\overline{D}_{i})$, $I_{h}^{(0)}:C(\overline{D}_{2})\rightarrow
V_{h}^{(0)}(\overline{D}_{2}),$ and the elliptic projection
$P_{1}:H^{1}(D_{1})\rightarrow V_{h}^{(1)}(D_{1})$ such that for $u\in
H^{1}(D_{1})$

$\int_{D_{1}}\left(k_{1}\frac{d(P_{1}u-u)}{dx}\frac{du_{h}}{dx}+\lambda(P_{1}u-u)u_{h}\right)dx=0\
\ \forall u_{h}\in V_{h}^{(1)}(D_{1}),$ (26)

where $\lambda>0$ is a constant; the error analysis for elliptic problems
suggests that a good choice is $\lambda=\underline{k_{1}}$. By virtue of (24)
it follows that there exists a constant $C$ independent of $h$ such that

$\left\|v-I_{h}^{(0)}v\right\|_{L^{2}(D_{2})}\leq
Ch\left\|v\right\|_{H^{1}(D_{2})};$ (27)

for $1\leq m\leq 2,\ 0\leq l\leq 1,$

$\left\|v-I_{h}^{(1)}v\right\|_{H^{l}(D_{i})}\leq
Ch^{m-l}\left\|v\right\|_{H^{m}(D_{i})};$ (28)

and from the well known error analysis for elliptic problems [2]

$\left\|u-P_{1}u\right\|_{H^{l}(D_{1})}\leq
Ch^{m-l}\left\|u\right\|_{H^{m}(D_{1})}.$ (29)

Likewise, for symmetric radial functions $v\in H_{r}^{q}(0,R)$, we define the
elliptic projector $P_{1}^{r}:H_{r}^{1}(0,R)\rightarrow V_{\Delta
r}^{(1)}[0,R]$ as the solution of the problem

$\int_{0}^{R}\left(k_{2}\frac{d(P_{1}^{r}w-w)}{dr}\frac{dw_{\Delta
r}}{dr}+\lambda(P_{1}^{r}w-w)w_{\Delta r}\right)r^{2}dr=0\ \ \forall w_{\Delta
r}\in V_{\Delta r}^{(1)}[0,R],$ (30)

with $\lambda>0$; as before, a good choice now is $\lambda=\underline{k_{2}}$.
By virtue of (25) it follows that there exists a constant $C$ independent of
$\Delta r$ such that

$\left\|w-P_{1}^{r}w\right\|_{L_{r}^{2}(0,R)}+\Delta
r\left\|\frac{d\left(w-P_{1}^{r}w\right)}{dr}\right\|_{L_{r}^{2}(0,R)}\leq
C\Delta r^{2}\left|w\right|_{H_{r}^{2}(0,R)}.$ (31)

$P_{1}^{r}$ can be extended to functions of $x$ and $r$ in an $L^{2}$-sense.
Thus, for $(x,r)\in D_{2}\times(0,R(\cdot))$ we define the extended projection
$P_{1}^{r}:H_{\ \ r}^{1,1}(D_{2}\times(0,R(\cdot)))\rightarrow
L^{2}(D_{2})\otimes V_{\Delta r}^{(1)}[0,R(\cdot)]$ as

$\int_{D_{2}}\int_{0}^{R(x)}\left(k_{2}\frac{\partial(P_{1}^{r}w-w)}{\partial
r}\frac{dw_{\Delta r}}{dr}+\lambda(P_{1}^{r}w-w)w_{\Delta
r}\right)r^{2}drdx=0\ \ \forall w_{\Delta r}\in L^{2}(D_{2})\otimes V_{\Delta
r}^{(1)}[0,R(\cdot)].$ (32)

Assuming that $w\in H_{\ \ r}^{1,1}(D_{2}\times(0,R(\cdot))$ is such that for
a.e. $x\in D_{2}$, $w(x,\cdot)\in H_{r}^{2}(0,R(x))$, then by virtue of (31)

$\left\|\left(w-P_{1}^{r}w\right)(x,\cdot)\right\|_{L_{r}^{2}(0,R(x))}+\Delta
r\left\|\frac{\left(\partial\left(w-P_{1}^{r}w\right)\right)(x,\cdot)}{\partial
r}\right\|_{L_{r}^{2}(0,R(x))}\leq C\Delta
r^{2}\left|w(x,\cdot)\right|_{H_{r}^{2}(0,R(x))}.$

Noting that $H_{\ \
r}^{0,q}(D_{2}\times(0,R(\cdot)))=L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))\cap
L^{2}(D_{2},H_{r}^{q}(0,R(\cdot)))$, then it readily follows that for $q=0,1$

$\left\|w-P_{1}^{r}w\right\|_{H_{\ \ r}^{0,q}(D_{2}\times(0,R(\cdot)))}\leq
C\Delta r^{2-q}\left\|w\right\|_{H_{\ \ r}^{0,2}(D_{2}\times(0,R(\cdot)))}.$
(33)

We shall also consider the $x$-Lagrange interpolant for functions that depend
on $x$ and $r$, $I_{0}^{x}:H_{\ \ r}^{1,1}(D_{2}\times(0,R(\cdot)))\rightarrow
V_{h}^{(0)}(\overline{D}_{2})\otimes H_{r}^{1}(0,R(\cdot))$. Thus, for
$v(x;r)\in H_{\ \ r}^{1,1}(D_{2}\times(0,R(\cdot)))$

$I_{0}^{x}v(x;r)=\sum_{l=1}^{M_{2}}v(x_{l};r)\chi_{l}(x)\ \textrm{with}\
r\in(0,R(x_{l})).$ (34)

Since $I_{0}^{x}$ can be viewed as an extended Lagrange interpolant
$I_{h}^{(0)}:C(\overline{D}_{2})\rightarrow V_{h}^{(0)}(\overline{D}_{2})$ ,
then based on (27) one can show that

$\left\|v-I_{0}^{x}v\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))}\leq
Ch\left\|v\right\|_{H^{1}(D_{2};L_{r}^{2}(0,R_{s}(\cdot)))}.$ (35)

###### Lemma 3

Let $I=(a,b),\ 0\leq a<b$, be a bounded interval with $\overline{I}=[a,b]$,
and let $f\in H^{1}(I)$. There exists an arbitrarily small number $\epsilon$
and a positive constant $C(\epsilon)$ such that

$\left\|f\right\|_{L^{\infty}(I)}\leq\epsilon\left\|\frac{df}{dx}\right\|_{L^{2}(I)}+C(\epsilon)\left\|f\right\|_{L^{2}(I)}.$
(36)

Proof. Since $H^{1}(I)\hookrightarrow C(\overline{I})$, then $f\in
C(\overline{I})$ and so does $f^{2}$, so for any $y,\ z\in I$, $y<z$, we have
that

$\begin{array}[]{r}f^{2}(z)-f^{2}(y)=\displaystyle\int_{y}^{z}\frac{df^{2}}{dx}dx=2\displaystyle\int_{y}^{z}f\frac{df}{dx}dx\\\
\\\
\leq\epsilon^{2}\displaystyle\int_{I}\left|\frac{df}{dx}\right|^{2}dx+\epsilon^{-2}\int_{I}f^{2}dx\text{.}\end{array}$

Since exists $x^{\ast}\in I$ such that $f^{2}(x^{\ast})=\min_{x\in
I}f^{2}(x)$, then letting $y=x^{\ast}$ it follows that

$f^{2}(y)\leq\frac{1}{\left|I\right|}\int_{I}f^{2}dx$

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^{\infty}-$ functions with respect to the $H_{r}^{1}(0,R)$-norm
and with the property that their with first derivative vanishes at $r=0$.

###### Lemma 4

If $v\in F^{1}(R),$ then for all $0<a<R,$

1) $v\in H^{1}(a,R).$

2) There exists an arbitrarily small number $\epsilon$ and a positive
(possibly large) constant $C(\epsilon)$, both depending on $a$, such that

$\left\|v\right\|_{L^{\infty}(a,R)}\leq\epsilon\left\|\frac{dv}{dr}\right\|_{L_{r}^{2}(0,R)}+C(\epsilon)\left\|v\right\|_{L_{r}^{2}(0,R)}.$
(37)

Proof. To prove 1) we note that for any $v\in F^{1}(R)$,

$\begin{array}[]{l}\left\|v\right\|_{H_{r}^{1}(0,R)}^{2}=\displaystyle\int_{0}^{R}v^{2}r^{2}dr+\displaystyle\int_{0}^{R}\left|\frac{dv}{dr}\right|^{2}r^{2}dr\\\
\\\
\geq\displaystyle\int_{a}^{R}v^{2}r^{2}dr+\displaystyle\int_{a}^{R}\left|\frac{dv}{dr}\right|^{2}r^{2}dr\\\
\\\ \geq
a^{2}\displaystyle\int_{a}^{R}v^{2}dr+\displaystyle\int_{a}^{R}\left|\frac{dv}{dr}\right|^{2}dr=a^{2}\left\|v\right\|_{H^{1}(a,R)}^{2}.\end{array}$
(38)

So, any sequence $\\{v_{n}\\}$ that converges with respect to the
$H_{r}^{1}(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).

###### Lemma 5

For each $(x,t)\in D_{2}\times[0,T_{\mathrm{end}}]$ we have the following
estimates.

$\begin{array}[]{l}\left|i_{0}-i_{0h}\right|\leq
C\left|u-u_{h}\right|+C\left|v_{s}-v_{hs}\right|,\\\ \\\ \left|\ln u-\ln
u_{h}\right|\leq C\left|u-u_{h}\right|,\\\ \\\
\left|\overline{U}(v_{s})-\overline{U}(v_{sh})\right|\leq
C\left|v_{s}-v_{sh}\right|.\end{array}$ (39)

Proof. Noting that the functions $x\rightarrow\sqrt{x},\
x\rightarrow\sqrt{1-x},\ x\rightarrow\ln x$ and $x\rightarrow\overline{U}(x)$
are smooth bounded and Lipschitz functions in any bounded interval $[a,b],\
0<a<b$, and that the composition and multiplication of bounded Lipschitz
functions results in a Lipschitz function, then the estimates follow. The
constant $C$ in (39) depends on the constants $P$ and $Q$ of (22).

###### Lemma 6

Let us consider $J$ and its approximate $J_{h}$, then for a.e.
$t\in[0,T_{end}]$ there exists a positive constant $C$ such that

$\begin{array}[]{r}\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}\leq
C\left\\{\left\|\phi_{2}(t)-\phi_{2h}(t)\right\|_{L^{2}(D_{2})}^{2}+\left\|\phi_{1}(t)-\phi_{1h}(t)\right\|_{L^{2}(D_{1})}^{2}\right.\\\
\\\
\left.+\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}\right\\}.\end{array}$
(40)

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
$(x,t)\in D_{2}\times[0,T_{\mathrm{end}}]$

$\begin{array}[]{l}\left|J-J_{h}\right|=\left|a_{2}i_{0}\sinh(\beta\eta)-a_{2}i_{0h}\sinh(\beta\eta_{h})\right|\\\
\\\
\leq\left|a_{2}i_{0}\left(\sinh(\beta\eta)-\sinh(\beta\eta_{h})\right)\right|+\left|a_{2}\left(i_{0}-i_{0h}\right)\sinh(\beta\eta_{h})\right|\\\
\\\ \leq
C\left|\sinh(\beta\eta)-\sinh(\beta\eta_{h})\right|+\left\|a_{2}\sinh(\beta\eta_{h})\right\|_{L^{\infty}(D_{2}\times[0,T_{end}])}\left|i_{0}-i_{0h}\right|\\\
\\\ \leq
C\left|\sinh(\beta\eta)-\sinh(\beta\eta_{h})\right|+C\left|i_{0}-i_{0h}\right|.\end{array}$

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\in(\eta,\eta_{h})$ such that

$\left|\sinh(\beta\eta)-\sinh(\beta\eta_{h})\right|\leq\beta\left|\cosh(z)\right|\left|\eta-\eta_{h}\right|\leq
C\left|\eta-\eta_{h}\right|,$

and resorting again to (22) and (23) it follows that

$\left\|z\right\|_{L^{\infty}(D_{2}\times[0,T_{\mathrm{end}}])}\leq\left\|\eta\right\|_{L^{\infty}(D_{2}\times[0,T_{\mathrm{end}}])}+\left\|\eta_{h}\right\|_{L^{\infty}(D_{2}\times[0,T_{\mathrm{end}}])}\leq
C$

Hence, applying Lemma 5 yields

$\left|J-J_{h}\right|\leq
C\left(\left|\phi_{2}-\phi_{2h}\right|+\left|\phi_{1}-\phi_{1h}\right|+\left|u-u_{h}\right|+\left|v_{s}-v_{sh}\right|\right)$

From this estimate it follows (40).

###  4.2 Error estimates for the potentials

To estimate the error for the potentials $\phi_{1}$ and $\phi_{2}$ is
convenient to introduce the spaces $V:=H^{1}(D_{1})\times
H^{1}(D_{2}):=\left\\{w=\left(w_{1},w_{2}\right):w_{1}\in H^{1}(D_{1}),\
w_{2}\in H^{1}(D_{2})\right\\}$ and $V_{h}\subset V$, where
$V_{h}:=V_{h}^{(1)}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$.
$V$ is a Hilbert space with norm

$\left\|w\right\|_{V}=\left(\left\|w_{1}\right\|_{H^{1}(D_{1})}^{2}+\left\|w_{2}\right\|_{H^{1}(D_{2})}^{2}\right)^{1/2},$

and seminorm

$\left|v\right|_{V}=\left(\left|w_{1}\right|_{H^{1}(D_{1})}^{2}+\left|w_{2}\right|_{H^{1}(D_{2})}^{2}\right)^{1/2}.$

Considering the bilinear forms $a_{1}:H^{1}(D_{1})\times
H^{1}(D_{1})\rightarrow\mathbb{R}$ and $a_{2}:H^{1}(D_{2})\times
H^{1}(D_{2})\rightarrow\mathbb{R}$,

$\left\\{\begin{array}[]{l}a_{1}(\phi_{1},\psi_{1})=\displaystyle\int_{D_{1}}\kappa(u)\frac{d\phi_{1}}{dx}\frac{d\psi_{1}}{dx}dx,\\\
\\\
a_{2}(\phi_{2},\psi_{2})=\displaystyle\int_{D_{2}}\sigma\frac{d\phi_{2}}{dx}\frac{d\psi_{2}}{dx}dx,\end{array}\right.$

we can define the bilinear form $a:V\times V\rightarrow\mathbb{R}$ as follows.
Let $\Phi$ and $\Psi$ $\in V$, $\Phi=(\phi_{1},\phi_{2})$ and
$\Psi=(\psi_{1},\psi_{2})$, then

$a(\Phi,\Psi)=a_{1}(\phi_{1},\psi_{1})+a_{2}(\phi_{2},\psi_{2}).$

Furthermore, concerning the right hand side terms of (12) and (13), we
introduce the operator $B:V\rightarrow V^{\ast}$, $V^{\ast}$ being the dual
for $V$, as

$\left\langle
B(\Phi),\Psi\right\rangle=\int_{D_{2}}J\left(\psi_{2}-\psi_{1}\right)dx.$

Hence, we can recast the equations (12) and (13) as follows. Find $\Phi\in
L^{2}(0,T_{end};W(D_{1})\times H^{1}(D_{2}))$ such that

$a(\Phi,\Psi)+\left\langle B(\Phi),\Psi\right\rangle=-\int_{D_{2}}g\psi_{2}dx\
\ \forall\Psi\in V.$ (41)

Likewise, the finite element solutions $\phi_{1h}$ and $\phi_{2h}$ that
satisfy (17) and (18) respectively, can be formulated as follows. For all
$t\in[0,T_{end}]$, find $\Phi_{h}\in W_{h}(\overline{D}_{1})\times
V_{h}^{(1)}(\overline{D}_{2})$ such that

$a_{h}(\Phi_{h},\Psi_{h})+\left\langle
B_{h}(\Phi_{h}),\Psi_{h}\right\rangle=-\int_{D_{2}}g\psi_{2h}dx\ \
\forall\Psi_{h}\in V_{h},$ (42)

where

$a_{h}(\Phi_{h},\Psi_{h})=a_{1h}(\phi_{1h},\psi_{1h})+a_{2}(\phi_{2h},\psi_{2h}),$
(43)

with

$a_{1h}(\phi_{1h},\psi_{1h})=\int_{D_{1}}\kappa(u_{h})\frac{d\phi_{1h}}{dx}\frac{d\psi_{1h}}{dx}dx,$

and

$\left\langle
B_{h}(\Phi_{h}),\Psi_{h}\right\rangle=\int_{D_{2}}J_{h}\left(\psi_{2h}-\psi_{1h}\right)dx.$
(44)

###### Remark 7

The existence and uniqueness of $\Phi$ 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 $\phi_{1}(x)\in 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.,

$\left\langle B(\Phi)-B(\widehat{\Phi}),\Phi-\widehat{\Phi}\right\rangle\geq
0,$ (45)

bounded and continuous in the sense that for all $\Psi\in V$ there exists a
constant $C$ such that

$\left\langle B(\Phi)-B(\widehat{\Phi}),\Psi\right\rangle\leq
C\left(\left\|\Phi-\widehat{\Phi}\right\|_{V}\right)\left(\left\|\psi_{1}\right\|_{L^{2}(D_{1})}+\left\|\psi_{2}\right\|_{L^{2}(D_{2})}\right).$
(46)

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

$\left\langle
B(\Phi)-B(\widehat{\Phi}),\Phi-\widehat{\Phi}\right\rangle=\int_{D_{2}}a_{2}i_{0}\left(\sinh(\beta\eta)-\sinh(\beta\widehat{\eta})\right)\left(\left(\phi_{2}-\widehat{\phi}_{2}\right)-\left(\phi_{1}-\widehat{\phi}_{1}\right)\right)dx,$

where $\eta=\phi_{2}-\phi_{1}-\alpha\ln u-\overline{U}$ and
$\widehat{\eta}=\widehat{\phi}_{2}-\widehat{\phi}_{1}-\alpha\ln
u-\overline{U}$. Since for all $x\in 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\in D_{2}$
$C(P,Q)\leq a_{2}i_{0}$, and by the mean value theorem
$\sinh(\beta\eta)-\sinh(\beta\widehat{\eta})\geq\beta(\eta-\widehat{\eta})$,
then one readily obtains

$\left\langle B(\Phi)-B(\widehat{\Phi}),\Phi-\widehat{\Phi}\right\rangle\geq
C\int_{D_{2}}\left(\left(\phi_{2}-\phi_{1}\right)-\left(\widehat{\phi}_{2}-\widehat{\phi}_{1}\right)\right)^{2}dx\geq
0.$ (47)

To prove that $B$ is bounded we notice that for all $\Phi\in V$

$\left\langle
B(\Phi),\Phi\right\rangle\leq\int_{D_{2}}\left|J\right|\times\left|\phi_{2}-\phi_{1}\right|dx=\int_{D_{2}}\left|a_{2}i_{0}(\sinh(\beta\eta)\right|\times\left|\phi_{2}-\phi_{1}\right|dx,$

but $\left|J\right|$ 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

$\left\langle B(\Phi),\Phi\right\rangle\leq C\left\|\Phi\right\|_{V},$

so $B$ is bounded. To prove that $B$ is continuous, we again notice that

$\left\langle
B(\Phi)-B(\widehat{\Phi}),\Psi\right\rangle\leq\int_{D_{2}}\left|a_{2}i_{0}(\sinh(\beta\eta)-\sinh(\beta\widehat{\eta}))\right|\times\left|\psi_{2}-\psi_{1}\right|dx,$

so, arguing as in the proof of Lemma 6 we have that there exists a positive
constant $C$ such that

$\left|a_{2}i_{0}(\sinh(\beta\eta)-\sinh(\beta\widehat{\eta}))\right|\leq
C\left|\eta-\widehat{\eta}\right|=C\left|\left(\phi_{2}-\widehat{\phi}_{2}\right)-\left(\phi_{1}-\widehat{\phi}_{1}\right)\right|.$

Substituting this estimate in the above inequality and making use of the
Cauchy-Schwarz inequality it follows that

$\begin{array}[]{c}\left\langle
B(\Phi)-B(\widehat{\Phi}),\Psi\right\rangle\leq
C\left(\left\|\phi_{2}-\widehat{\phi}_{2}\right\|_{L^{2}(D_{2})}+\left\|\phi_{1}-\widehat{\phi}_{1}\right\|_{L^{2}(D_{1})}\right)\left(\left\|\psi_{2}\right\|_{L^{2}(D_{2})}+\left\|\psi_{1}\right\|_{L^{2}(D1)}\right)\\\
\\\ \leq
C\left\|\Phi-\widehat{\Phi}\right\|_{V}\left(\left\|\psi_{2}\right\|_{L^{2}(D_{2})}+\left\|\psi_{1}\right\|_{L^{2}(D1)}\right).\end{array}$

###### 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\in[0,T_{\mathrm{end}}]$, let the solution of (41),
$\Phi(t)=\left(\phi_{1}(t),\phi_{2}(t)\right)$, be in $H^{2}(D_{1})\times
H^{2}(D_{2})$. There exists a constant $C$ independent of $h$ such that

$\begin{array}[]{c}\left\|\Phi(t)-\Phi_{h}(t)\right\|_{V}^{2}\leq
C\left(h^{2}(\left\|\phi_{1}(t)\right\|_{H^{2}(D_{1})}^{2}+\left\|\phi_{2}(t)\right\|_{H^{2}(D_{2})}^{2})\right.\\\
\\\
\left.+\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}\right).\end{array}$
(48)

Proof. Setting $\Psi=\Psi_{h}$ in (41) and subtracting (42) yields

$a(\Phi,\Psi_{h})-a_{h}(\Phi_{h},\Psi_{h})+\left\langle
B(\Phi)-B\left(\Phi_{h}\right),\Psi_{h}\right\rangle=-\left\langle
B(\Phi_{h})-B_{h}\left(\Phi_{h}\right),\Psi_{h}\right\rangle.$

Noting that

$a(\Phi,\Psi_{h})-a_{h}(\Phi_{h},\Psi_{h})=a_{h}(\Phi-\Phi_{h},\Psi_{h})+\int_{D_{1}}\left(\kappa(u)-\kappa(u_{h})\right)\frac{\partial\phi_{1}}{\partial
x}\frac{d\psi_{1h}}{dx}dx,$

it follows that

$\begin{array}[]{r}a_{h}(\Phi-\Phi_{h},\Psi_{h})+\left\langle
B(\Phi)-B\left(\Phi_{h}\right),\Psi_{h}\right\rangle=\left\langle
B_{h}(\Phi_{h})-B\left(\Phi_{h}\right),\Psi_{h}\right\rangle\\\ \\\
-\displaystyle\int_{D_{1}}\left(\kappa(u)-\kappa(u_{h})\right)\frac{\partial\phi_{1}}{\partial
x}\frac{d\psi_{1h}}{dx}dx.\end{array}$ (49)

To estimate the terms of this expression we choose
$\Psi_{h}=I_{h}^{(1)}\Phi-\Phi_{h}=$
$\left(I_{h}^{(1)}\phi_{1}-\phi_{1h},I_{h}^{(1)}\phi_{2}-\phi_{2h}\right)$,
$I_{h}^{(1)}$ being the Lagrange interpolant on
$V_{h}=V_{h}^{(1)}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$,
this means that $I_{h}^{(1)}\phi_{1}\in V_{h}^{(1)}(\overline{D}_{1})$ and
$I_{h}^{(1)}\phi_{2}\in V_{h}^{(1)}(\overline{D}_{2})$. For convenience, we
shall split the expression for $\Psi_{h}$ as

$\Psi_{h}=\left(\Phi-\Phi_{h}\right)+\left(I_{h}^{(1)}\Phi-\Phi\right).$

Replacing this expression for $\Psi_{h}$ in (49) we have that

$\left\\{\begin{array}[]{l}a_{h}(\Phi-\Phi_{h},\Phi-\Phi_{h})+\left\langle
B(\Phi)-B\left(\Phi_{h}\right),\Phi-\Phi_{h}\right\rangle\\\ \\\
=a_{h}(\Phi-\Phi_{h},\Phi-I_{h}^{(1)}\Phi)+\left\langle
B(\Phi)-B\left(\Phi_{h}\right),\Phi-I_{h}^{(1)}\Phi\right\rangle\\\ \\\
+\left\langle
B_{h}(\Phi_{h})-B\left(\Phi_{h}\right),\left(\Phi-\Phi_{h}\right)+\left(I_{h}^{(1)}\Phi-\Phi\right)\right\rangle\\\
\\\
+\displaystyle\int_{D_{1}}\left(\kappa(u)-\kappa(u_{h})\right)\frac{\partial\phi_{1}}{\partial
x}\frac{\partial\left(\left(\phi_{1}-\phi_{1h}\right)+\left(I_{h}^{(1)}\phi_{1}-\phi_{1h}\right)\right)}{\partial
x}dx.\end{array}\right.$ (50)

We bound the terms of (50). We start by showing that there exists a positive
constant $\alpha$ such that the term on the left hand side satisfies

$a_{h}(\Phi-\Phi_{h},\Phi-\Phi_{h})+\left\langle
B(\Phi)-B\left(\Phi_{h}\right),\Phi-\Phi_{h}\right\rangle\geq\alpha\left\|\Phi-\Phi_{h}\right\|_{V}^{2}.$
(51)

To do so we note that by virtue of (47)

$\begin{array}[]{l}a_{h}(\Phi-\Phi_{h},\Phi-\Phi_{h})+\left\langle
B(\Phi)-B\left(\Phi_{h}\right),\Phi-\Phi_{h}\right\rangle\\\ \\\ \geq
a_{1h}(\phi_{1}-\phi_{1h},\phi_{1}-\phi_{1h})+a_{2}(\phi_{2}-\phi_{2h},\phi_{2}-\phi_{2h})\\\
\\\
+C\displaystyle\int_{D_{2}}\left(\left(\phi_{2}-\phi_{1}\right)-\left(\phi_{2h}-\phi_{1h}\right)\right)^{2}dx.\end{array}$
(52)

Since $\phi_{1}-\phi_{1h}\in W(D_{1})$, we can use A2 and Poincaré-Wirtinger
inequality to bound $a_{1}(\phi_{1}-\phi_{1h},\phi_{1}-\phi_{1h})$ from below
as

$a_{1h}(\phi_{1}-\phi_{1h},\phi_{1}-\phi_{1h})\geq
c_{1}\left\|\phi_{1}-\phi_{1h}\right\|_{H^{1}(D_{1})}^{2},$

where the constant $c_{1}=\kappa_{0}\left(1+C_{P}\right)^{-1}$, $C_{P}$ being
the constant of the Poincaré-Wirtinger inequality; using again A2, we bound
the term $a_{2}(\phi_{2}-\phi_{2h},\phi_{2}-\phi_{2h})$ as

$a_{2}(\phi_{2}-\phi_{2h},\phi_{2}-\phi_{2h})\geq\sigma_{0}\left|\phi_{2}-\phi_{2h}\right|_{H^{1}(D_{2})}^{2}.$

Applying Young inequality we find that there exists a constant
$\gamma\in\left(0,1\right)$ such that

$\begin{array}[]{l}C\displaystyle\int_{D_{2}}\left(\left(\phi_{2}-\phi_{1}\right)-\left(\phi_{2h}-\phi_{1h}\right)\right)^{2}dx.\\\
\\\ \geq
C\left[(1-\frac{4}{\gamma})\left\|\phi_{1}-\phi_{1h}\right\|_{L^{2}(D_{2})}^{2}+(1-\gamma)\left\|\phi_{2}-\phi_{2h}\right\|_{L^{2}(D_{2})}^{2}\right]\\\
\\\ \geq
C\left[(1-\frac{4}{\gamma})\left\|\phi_{1}-\phi_{1h}\right\|_{L^{2}(D_{1})}^{2}+(1-\gamma)\left\|\phi_{2}-\phi_{2h}\right\|_{L^{2}(D_{2})}^{2}\right].\end{array}$

Now, we can choose the constants $C$ and $\gamma$ such that
$c_{2}=c_{1}+C(1-\frac{4}{\gamma})>0$, and substitute these bounds in (52) to
obtain the inequality (51), where
$\alpha=\min(c_{1},c_{2},\sigma_{0},(1-\gamma)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 $\epsilon_{1}$ and a
constant $C(\epsilon_{1})$ such that

$\begin{array}[]{r}a_{h}(\Phi-\Phi_{h},\Phi-I_{h}^{(1)}\Phi)\leq
C\left\|\Phi-\Phi_{h}\right\|_{V}\left\|\Phi-I_{h}^{(1)}\Phi\right\|_{V}\\\
\\\
\leq\epsilon_{1}\left\|\Phi-\Phi_{h}\right\|_{V}^{2}+C(\epsilon_{1})\left\|\Phi-
I_{h}^{(1)}\Phi\right\|_{V}^{2}.\end{array}$ (53)

To bound $\left\langle B(\Phi)-B\left(\Phi_{h}\right),\Phi-
I_{h}^{(1)}\Phi\right\rangle$ we note that

$\left\langle B(\Phi)-B\left(\Phi_{h}\right),\Phi-
I_{h}^{(1)}\Phi\right\rangle\leq\int_{D_{2}}\left|a_{2}i_{0}\right|\left|\int_{\beta\eta_{h}}^{\beta\eta}\cosh\xi
d\xi\right|\left|\Phi-I_{h}^{(1)}\Phi\right|dx.$

By virtue of (22) and (23) and the mean value theorem for the integral

$\begin{array}[]{l}\left\langle B(\Phi)-B\left(\Phi_{h}\right),\Phi-
I_{h}^{(1)}\Phi\right\rangle\leq
C\displaystyle\int_{D_{2}}\left|\eta-\eta_{h}\right|\left|\Phi-
I_{h}^{(1)}\Phi\right|dx\\\ \\\ \leq
C\displaystyle\int_{D_{2}}\left(\left|\phi_{2}-\phi_{2h}\right|+\left|\phi_{1}-\phi_{1h}\right|+\left|u-u_{h}\right|+\left|v_{s}-v_{sh}\right|\right)\left|\Phi-
I_{h}^{(1)}\Phi\right|dx.\end{array}$

Applying Young inequality yields

$\begin{array}[]{c}\left\langle B(\Phi)-B\left(\Phi_{h}\right),\Phi-
I_{h}^{(1)}\Phi\right\rangle\leq\epsilon_{2}\left(\left\|\Phi-\Phi_{h}\right\|_{V}^{2}+\left\|u-u_{h}\right\|_{L^{2}(D_{1})}^{2}\right.\\\
\\\
\left.+\left\|v_{s}-v_{sh}\right\|_{L^{2}(D_{2})}^{2}\right)+C(\epsilon_{2})\left\|\Phi-
I_{h}^{(1)}\Phi\right\|_{V}^{2},\end{array}$ (54)

where $\epsilon_{2}$ is a small positive number and $C(\epsilon_{2})$ is a
constant. Next, noting that by virtue of (22) and (23)
$\left|\sinh(\beta\eta_{h})\right|$ is bounded in $D_{2}$, then

$\left\langle
B_{h}(\Phi_{h})-B\left(\Phi_{h}\right),\left(\Phi-\Phi_{h}\right)+\left(I_{h}^{(1)}\Phi-\Phi\right)\right\rangle\leq
C\int_{D_{2}}a_{2}\left|i_{0}-i_{0h}\right|\left|\left(\Phi-\Phi_{h}\right)+\left(I_{h}^{(1)}\Phi-\Phi\right)\right|dx.$

Again, using Lemma 5 and Young inequality we obtain that there exist a small
number $\epsilon_{3}$ and a constant $C(\epsilon_{3})$ such that

$\begin{array}[]{l}\left\langle
B_{h}(\Phi_{h})-B\left(\Phi_{h}\right),\left(\Phi-\Phi_{h}\right)+\left(I_{h}^{(1)}\Phi-\Phi\right)\right\rangle\leq\epsilon_{3}\left\|\Phi-\Phi_{h}\right\|_{V}^{2}\\\
\\\
+\epsilon_{3}\left(\left\|u-u_{h}\right\|_{L^{2}(D_{1})}^{2}+\left\|v_{s}-v_{sh}\right\|_{L^{2}(D_{2})}^{2}\right)+C(\epsilon_{3})\left\|\Phi-
I_{h}^{(1)}\Phi\right\|_{V}^{2}.\end{array}$ (55)

To bound the last term on the right hand side of (50) we note that
$\left\|\frac{\partial\phi_{1}}{\partial x}\right\|_{L^{\infty}(D_{1})}$ is
bounded and by virtue of assumption A4, and by the mean value theorem,
$\left|\left(\kappa(u)-\kappa(u_{h})\right)\right|\leq C\left|u-u_{h}\right|$,
then it follows that

$\begin{array}[]{l}\displaystyle\int_{D_{1}}\left(\kappa(u)-\kappa(u_{h})\right)\frac{\partial\phi_{1}}{\partial
x}\frac{\partial\left(\left(\phi_{1}-\phi_{1h}\right)+\left(I_{h}^{(1)}\phi_{1}-\phi_{1h}\right)\right)}{\partial
x}dx\\\ \\\ \leq
C\displaystyle\int_{D_{1}}\left|u-u_{h}\right|\left|\frac{\partial\left(\left(\phi_{1}-\phi_{1h}\right)+\left(I_{h}^{(1)}\phi_{1}-\phi_{1h}\right)\right)}{\partial
x}\right|dx\\\ \\\
\leq\epsilon_{4}\left\|\Phi-\Phi_{h}\right\|_{V}^{2}+C(\epsilon_{4})\left(\left\|u-u_{h}\right\|_{L^{2}(D_{1})}^{2}+\left\|\Phi-
I_{h}^{(1)}\Phi\right\|_{V}^{2}\right).\end{array}$ (56)

Letting$\ \epsilon_{1}+\cdots+\epsilon_{4}=\alpha/2$ and noting that, see
(34),

$\begin{array}[]{c}\left\|\Phi-
I_{h}^{(1)}\Phi\right\|_{V}^{2}=\left\|\phi_{1}-I_{h}^{(1)}\phi_{1}\right\|_{H^{1}(D_{1})}^{2}+\left\|\phi_{2}-I_{h}^{(1)}\phi_{2}\right\|_{H^{1}(D_{2})}^{2}\\\
\\\ \leq
Ch^{2}\left(\left\|\phi_{1}\right\|_{H^{2}(D_{1})}^{2}+\left\|\phi_{2}\right\|_{H^{2}(D_{2})}^{2}\right).\end{array}$

the estimate (48) follows from (50)-(56).

###  4.3 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\Delta 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

$u-u_{h}=(u-P_{1}u)+(P_{1}u-u_{h})\equiv\rho_{u}+\theta_{u},$ (57)

where $P_{1}$ is the elliptic projector defined in (26), and note that
$\theta_{u}(x,t)\in V_{h}^{(1)}(\overline{D}_{1})$. To carry out a
decomposition of this kind for $v(x;r,t)-v_{h\Delta r}(x;r,t)$, at first we
can try using the extended elliptic projector $P_{1}^{r}:H_{\ \
r}^{1,1}(D_{2}\times(0,R_{s}(\cdot)))\rightarrow L^{2}(D_{2})\otimes V_{\Delta
r}^{(1)}[0,R_{s}(\cdot)]$ defined in (32) and assume that for all $t$,
$v(x;r,t)\in H_{\ \ r}^{1,1}(D_{2}\times(0,R_{s}(\cdot)))$, then we find that
for a. e. $x\in D_{2}$

$P_{1}^{r}v(x;r,t)=\sum_{j=1}^{M}P_{1}^{r}v(x;r_{j},t)\alpha_{j}(r),$

here, $M$ denotes the number of mesh points in $[0,R_{s}(\cdot)]$,
$\left\\{\alpha_{j}(r)\right\\}_{j=1}^{M}$ the nodal basis of the linear
finite element space $V_{\Delta r}^{(1)}[0,R_{s}(\cdot)]\subset
H_{r}^{1}(0,R_{s}(\cdot))$ and the function $P_{1}^{r}v(x;r_{j},t)\in
L^{2}(D_{2})$; so, in general, $P_{1}^{r}v(x;r,t)$ is not in $V_{h\Delta
r}(\overline{D}_{3})$ and, consequently, it does not make sense to use
$P_{1}^{r}v(x;r,t)-v_{h\Delta r}(x;r,t)$ for such type of decomposition;
however, recalling the interpolant $I_{0}^{x}:H_{\ \
r}^{1,1}(D_{2}\times(0,R(\cdot)))\rightarrow
V_{h}^{(0)}(\overline{D}_{2})\otimes H_{r}^{1}(0,R(\cdot))$ defined in (34)
and further assuming that $P_{1}^{r}v(x;r_{j},t)\in H^{1}(D_{2})$, then it
follows that $P_{1}^{r}v(x;r,t)\in H_{\ \ r}^{1,1}(D_{2}\times(0,R(\cdot)))$
and, therefore, we can define $I_{0}^{x}P_{1}^{r}v(x;r,t)$ as

$I_{0}^{x}P_{1}^{r}v(x;r,t)=\sum_{l=1}^{M_{2}}\sum_{j=1}^{M^{(l)}}P_{1}^{r}v(x_{l};r_{j},t)\alpha_{j}^{(l)}(r)\chi_{l}(x);$

this expression implies that $I_{0}^{x}P_{1}^{r}v(x;r,t)\in V_{h\Delta
r}(\overline{D}_{3})$, so it makes sense to set

$\theta_{v}(x;r,t)=I_{0}^{x}P_{1}^{r}v(x;r,t)-v_{h\Delta r}(x;r,t).$

Now, using again the extended $P_{1}^{r}$ elliptic projector we define

$\rho_{v}(x;r,t)=v(x;r,t)-P_{1}^{r}v(x;r,t)\in H_{\ \
r}^{1,1}(D_{2}\times(0,R(\cdot))),$

and consequently

$I_{0}^{x}\rho_{v}(x;r,t)=I_{0}^{x}v(x;r,t)-I_{0}^{x}P_{1}^{r}v(x;r,t).$

Then, from all these considerations we can write that

$v(x;r,t)-v_{h\Delta
r}(x;r,t)=(v-I_{0}^{x}v)(x;r,t)+I_{0}^{x}\rho_{v}(x;r,t)+\theta_{v}(x;r,t).$
(58)

From (57) and (58) it follows that

$\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}\leq\left\|\rho_{u}(t)\right\|_{L^{2}(D_{1})}+\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})},$

and

$\begin{array}[]{r}\left\|v(t)-v_{h\Delta
r}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))}\leq\left\|v(t)-I_{0}^{x}v(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))}\\\
\\\
+\left\|I_{0}^{x}\rho_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))}.\end{array}$

The estimates for $\rho_{u}$ and $\rho_{v}$ are given in (29) and (33)
respectively, i.e.,

$\left\\{\begin{array}[]{l}\left\|\rho_{u}(t)\right\|_{L^{2}(D_{1})}\leq
Ch^{2}\left\|u(t)\right\|_{H^{2}(D_{1})},\\\ \\\
\left\|\rho_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R(\cdot)))}\leq C\Delta
r^{2}\left\|v(t)\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R(\cdot)))},\end{array}\right.$
(59)

then, it remains to calculate the estimates for $\theta_{u}$ and $\theta_{v}$;
but before going into the details of such calculations, we present new
estimates for $\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}$ and
$\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}$, which depend on
$\theta_{u}$ and $\theta_{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 $\epsilon$ and constants $C$ and
$C(\epsilon)$ independent of $h$ and $\Delta r$ such that

$\begin{array}[]{r}\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}\leq
Ch^{2}\left|v_{s}(t)\right|_{H^{1}(D_{2})}^{2}+C(\epsilon)\Delta
r^{2}\left\|v(t)\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}\\\ \\\
+C(\epsilon)\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\displaystyle\left\|\frac{\partial\theta_{v}(t)}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2},\end{array}$ (60)

and

$\begin{array}[]{c}\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}\leq
Ch^{2}\left(\left\|\phi_{1}(t)\right\|_{H^{2}(D_{1})}^{2}+\left\|\phi_{2}(t)\right\|_{H^{2}(D_{2})}^{2}+h^{2}\left\|u(t)\right\|_{H^{2}(D_{1})}^{2}\right)\\\
\\\
+C\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}+C\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}.\end{array}$
(61)

Proof. To calculate the estimate (60) we set
$v_{s}-v_{sh}=\left(v_{s}-I_{h}^{0}v_{s}\right)+\left(I_{h}^{0}v_{s}-v_{sh}\right)$,
so using (27) it follows that

$\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}=Ch^{2}\left|v_{s}(t)\right|_{H^{1}(D_{2})}^{2}+2\left\|I_{h}^{0}v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}.$

Since $I_{h}^{0}v_{s}-v_{sh}=\left(I_{h}^{0}v-v_{h\Delta
r}\right)(x;R_{s}(x),t)=\left(I_{0}^{x}\rho_{v}+\theta_{v}\right)(x;R_{s}(x),t)$
(recalling the definition of $\theta_{v}$) then by virtue of the definition of
the $L^{2}$-norm for functions of $V_{h}^{0}(\overline{D}_{2})$ presented in
Section 3, we have that

$\left\|I_{h}^{0}v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}\leq
2\sum_{l=1}^{M_{2}}\widehat{h}_{l}\left(\rho_{v}^{2}(x_{l};R_{s}(x_{l}),t)+\theta_{v}^{2}(x_{l};R_{s}(x_{l}),t)\right).$

To estimate $\rho_{v}^{2}(x_{l};R_{s}(x_{l}),t)$ and
$\theta_{v}^{2}(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
$\rho_{v}^{2}(x_{l};R_{s}(x_{l}),t)\leq\left\|\rho_{v}(x_{l};r,t)\right\|_{L^{\infty}(a,R(x_{l}))}^{2}$
and
$\theta_{v}^{2}(x_{l};R_{s}(x_{l}),t)\leq\left\|\theta_{v}(x_{l};r,t)\right\|_{L^{\infty}(a,R(x_{l}))}^{2}$,
thus by virtue of (37) it follows that there are a real number $\epsilon$ and
positive constant $C(\epsilon)>\epsilon$ such that

$\begin{array}[]{r}\displaystyle\sum_{l=1}^{M_{2}}\widehat{h}_{l}\rho_{v}^{2}(x_{l};R_{s}(x_{l}),t)\leq\displaystyle\epsilon\sum_{l=1}^{M_{2}}\widehat{h}_{l}\int_{0}^{R_{s}(x_{l})}r^{2}\left|\frac{\partial\rho_{v}(x_{l};r,t)}{\partial
r}\right|^{2}dr\\\ \\\ +\displaystyle
C(\epsilon)\sum_{l=1}^{M_{2}}\widehat{h}_{l}\int_{0}^{R_{s}(x_{l})}r^{2}\left|\rho_{v}(x_{l};r,t)\right|^{2}dr\\\
\\\ \leq
C(\epsilon)\left\|I_{0}^{x}\rho_{v}(t)\right\|_{L^{2}(D_{2},H_{r}^{1}(0,R_{s}(\cdot)))}^{2}\\\
\\\ \leq
C(\epsilon)\left\|\rho_{v}(t)\right\|_{L^{2}(D_{2},H_{r}^{1}(0,R_{s}(x)))}\\\
\\\ \leq C(\epsilon)\Delta
r^{2}\left\|v(t)\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2},\end{array}$

because by approximation theory
$\left\|I_{0}^{x}\rho_{v}(t)\right\|_{L^{2}(D_{2},H_{r}^{1}(0,R_{s}(x)))}\leq
C\left\|\rho_{v}(t)\right\|_{L^{2}(D_{2},H_{r}^{1}(0,R_{s}(x)))}$, this latter
term being estimated according to (33). Similarly,

$\begin{array}[]{r}\displaystyle\sum_{l=1}^{M_{2}}\widehat{h}_{l}\theta_{v}^{2}(x_{l};R_{s}(x_{l}),t)\leq\epsilon\displaystyle\sum_{l=1}^{M_{2}}\widehat{h}_{l}\int_{0}^{R_{s}(x_{l})}r^{2}\left|\frac{\partial\theta_{v}(x_{l};r,t)}{\partial
r}\right|^{2}dr\\\ \\\
+C(\epsilon)\displaystyle\sum_{l=1}^{M_{2}}\widetilde{h}_{l}\int_{0}^{R_{s}(x_{l})}r^{2}\left|\theta_{v}(x_{l};r,t)\right|^{2}dr\\\
\\\
\leq\epsilon\left|\theta_{v}(t)\right|_{L^{2}(D_{2},H_{r}^{1}(0,R_{s}(\cdot)))}^{2}+C(\epsilon)\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.\end{array}$
(62)

So, putting these bounds together the result (60) follows. To calculate the
estimate (61) we notice that by virtue of (40) and Theorem 10

$\begin{array}[]{r}\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}\leq
C\left(h^{2}\left(\left\|\phi_{1}(t)\right\|_{H^{2}(D_{1})}^{2}+\left\|\phi_{2}(t)\right\|_{H^{2}(D_{2})}^{2}\right)\right.\\\
\\\
\left.+\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}^{2}\right).\end{array}$

Since $u-u_{h}=\rho_{u}+\theta_{u}$, then taking into account (29)

$\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}\leq
Ch^{4}\left\|u(t)\right\|_{H^{2}(D_{1})}^{2}+2\left\|\theta_{u}\right\|_{L^{2}(D_{1})}^{2},$

so the result (61) follows.

Next, we calculate an estimate for $\theta_{v}$. To this end, we obtain, based
on equation (11), the integral equation for $I_{0}^{x}v(x;r,t)$ that will be
used for that purpose. Thus, for each one of the mesh points
$\\{x_{l}\\}_{l=1}^{M_{2}}$ of $D_{2h}$ the equation (11) reads

$\int_{0}^{R_{s}(x_{l})}\frac{\partial{v}^{(l)}}{\partial{t}}w^{(l)}r^{2}dr+\int_{0}^{R_{s}(x_{l})}k_{2}\frac{\partial
v^{(l)}}{\partial r}\frac{\partial w^{(l)}}{\partial
r}r^{2}dr=-\left(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}Jw^{(l)}(R_{s}(x))\right)|_{x=x_{l}},$
(63)

where $v^{(l)}:=v(x_{l};r,t)\in H_{r}^{1}(0,R_{s}(x_{l})$ and
$w^{(l)}:=w(x_{l};r)\in H_{r}^{1}(0,R_{s}(x_{l}))$. Using the nodal basis
functions $\left\\{\chi_{l}(x)\right\\}_{l=1}^{M_{2}}$ of the finite element
space $V_{h}^{(0)}(\overline{D}_{2})$ we can write

$I_{0}^{x}{v(x;r,t)=}\sum_{l}^{M_{2}}v^{(l)}(r,t)\chi_{l}(x)\mathrm{\ and\ \
}w(x;r)=\sum_{l=1}^{M_{2}}w^{(l)}(r)\chi_{l}(x).$

Now, noting that

$\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(I_{0}^{x}{v}\right)wr^{2}drdx=\sum_{l=1}^{M_{2}}\int_{D_{2}}\chi_{l}^{2}(x)\left(\int_{0}^{R_{s}(x)}v^{(l)}(r,t)w^{(l)}(r)r^{2}dr\right)dx,$

and for $x\in\widehat{e}_{l}$, $R_{s}(x)=R_{s}(x_{l})$, then it follows that

$\sum_{l=1}^{M_{2}}\int_{D_{2}}\chi_{l}^{2}(x)\left(\int_{0}^{R_{s}(x)}v^{(l)}(r,t)w^{(l)}(r)r^{2}dr\right)dx=\sum_{l=1}^{M_{2}}\widehat{h}_{l}\int_{0}^{R_{s}(x_{l})}v^{(l)}(r,t)w^{(l)}(r)r^{2}dr.$

Hence, (63) becomes

$\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(\frac{\partial
I_{0}^{x}{v}}{\partial{t}}w+k_{2}\frac{\partial I_{0}^{x}v}{\partial
r}\frac{\partial w}{\partial
r}\right)r^{2}drdx=-\int_{D_{2}}I_{0}^{x}(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}Jw(x;R_{s}(x)))dx.$
(64)

We proceed to formulate the equation for $\theta_{v}$. From (58) it follows
that $v_{h\Delta r}=I_{0}^{x}v-(I_{0}^{x}\rho_{v}+\theta_{v})$, then replacing
this expression for $v_{h\Delta r}$ in (16) and using (32) and (64) it follows
that for all $w_{h\Delta r}\in V_{h\Delta r}(\overline{D}_{3})$,

$\begin{array}[]{r}\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(\frac{\partial\theta_{v}}{\partial{t}}w_{h\Delta
r}+k_{2}\frac{\partial\theta_{v}}{\partial r}\frac{\partial w_{h\Delta
r}}{\partial
r}\right)r^{2}drdx=\lambda\int_{D_{2}}\int_{0}^{R_{s}(x)}I_{0}^{x}\rho_{v}w_{h\Delta
r}r^{2}drdx\\\ \\\ -\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\frac{\partial
I_{0}^{x}\rho_{v}}{\partial{t}}w_{h\Delta r}r^{2}drdx\\\ \\\
-\displaystyle\int_{D_{2}}I_{h}^{0}(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}(J(x)-J_{h}(x))w_{sh}(x))dx\\\
\\\
+\displaystyle\int_{D_{2}}\left(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}J_{h}(x)-I_{h}^{0}(R_{s}^{2}(x)a_{2}^{-1}F^{-1}(x)J_{h}(x)\right))w_{sh}(x)dx,\end{array}$
(65)

where we have made use of the following properties of the interpolant
$I_{0}^{x}$: (i) for $w_{h\Delta r}(x,r)\in V_{h\Delta r}(\overline{D}_{3})$,
$w_{h\Delta r}(x,r)=I_{0}^{x}w_{h\Delta r}(x,r)$, and (ii) when $r=R_{s}(x)$,
we can define the function $w_{sh}(x)=w_{h\Delta r}(x,R_{s}(x))$ such that
$w_{sh}(x)=I_{0}^{x}w_{h\Delta r}(x,R_{s}(x))=I_{h}^{0}w_{sh}(x)$. Setting
$w_{h\Delta r}=\theta_{v}$ yields

$\begin{array}[]{l}\displaystyle\frac{1}{2}\frac{d}{dt}\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\displaystyle\underline{k_{2}}\left\|\frac{\partial\theta_{v}(t)}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\\\ \\\
\leq\lambda\left\|I_{0}^{x}\rho_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\\\
\\\ +\displaystyle\left\|\frac{\partial I_{0}^{x}\rho_{v}(t)}{\partial
t}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\\\
\\\
+C\left\|I_{h}^{(0)}\left((J-J_{h})\theta_{vs}(t)\right)\right\|_{L^{1}(D_{2})}+C\left\|(J_{h}-I_{(0)}^{0}J_{h})\theta_{vs}(t)\right\|_{L^{1}(D_{2})}\equiv\sum_{i=1}^{4}R_{i},\end{array}$
(66)

where, $\theta_{vs}(t)=\theta_{v}(x;R_{s}(x),t)$ is the value of $\theta_{v}$
on the surface of the sphere of radius $R_{x}(x)$ associated with the point
$\\{x\\}$ of $D_{2}$.

###### Lemma 12

There exists a constant $C$ independent of $h$ and $\Delta r$, but depending
on $\underline{k_{1}}$ and $\underline{k_{2}}$, such that

$\begin{array}[]{c}\displaystyle\frac{d}{dt}\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\displaystyle\underline{k_{2}}\left\|\frac{\partial\theta_{v}(t)}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\\\ \\\ \leq
Ch^{2}\left(\left\|\phi_{1}(t)\right\|_{H^{2}(D_{1})}^{2}+\left\|\phi_{2}(t)\right\|_{H^{2}(D_{2})}^{2}\right)\\\
\\\
+Ch^{2}\left(\left|v_{s}(t)\right|_{H^{1}(D_{2})}^{2}+h^{2}\left\|u(t)\right\|_{H^{2}(D_{1})}^{2}+\displaystyle\left\|\frac{\partial
J}{\partial x}\right\|_{L^{2}(D_{2})}^{2}\right)\\\ \\\ +C\Delta
r^{2}\left(\left\|v(t)\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\Delta
r^{2}\displaystyle\left\|\frac{\partial v(t)}{\partial
t}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right)\\\ \\\
+C\left(\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\right).\end{array}$
(67)

Proof. We bound the terms $R_{1},\ldots,R_{4}$ on the right hand side of (66).
Noting that
$\left\|I_{0}^{x}\rho_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\leq
C\left\|\rho_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}$, then by
virtue of Young inequality and (33) it follows that

$R_{1}\leq C\Delta
r^{4}\left\|v(t)\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}+C\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.$
(68)

To bound $R_{2}$ we notice that

$\begin{array}[]{r}\displaystyle\left\|\frac{\partial
I_{0}^{x}\rho_{v}(t)}{\partial
t}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\leq
C\displaystyle\left\|\frac{\partial\rho_{v}(t)}{\partial
t}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\\\ \\\
=C\displaystyle\left\|(I-P_{1}^{r})\frac{\partial v(t)}{\partial
t}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\\\ \\\ \leq C\Delta
r^{2}\displaystyle\left\|\frac{\partial v(t)}{\partial
t}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}.\end{array}$

Hence, by Young inequality it follows that

$R_{2}\leq C\Delta r^{4}\left\|\frac{\partial v(t)}{\partial
t}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}+C\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.$
(69)

We bound the term $R_{3}$. Thus, we have that

$\begin{array}[]{l}C\left\|I_{h}^{0}\left((J-J_{h})\theta_{vs}(t)\right)\right\|_{L^{1}(D_{2})}=C\displaystyle\sum_{l=1}^{M_{2}}\widehat{h}_{l}\left|\left(J(x_{l})-J_{h}(x_{l}\right))\theta_{v}(x_{l};R_{s}(x_{l}),t)\right|\\\
\\\
\leq\displaystyle\frac{C^{2}}{2}\sum_{l=1}^{M_{2}}\widehat{h}_{l}\left(J(x_{l})-J_{h}(x_{l}\right))^{2}+\displaystyle\frac{1}{2}\sum_{l=1}^{M_{2}}\widehat{h}_{l}\theta_{v}^{2}(x_{l};R_{s}(x_{l}),t)\\\
\\\
\leq\displaystyle\frac{C^{2}}{2}\left\|I_{h}^{0}(J-J_{h})\right\|_{L^{2}(D_{2})}^{2}+\displaystyle\frac{1}{2}\sum_{l=1}^{M_{2}}\widehat{h}_{l}\theta_{v}^{2}(x_{l};R_{s}(x_{l}),t).\end{array}$

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
$\left\|I_{h}^{0}(J-J_{h})\right\|_{L^{2}(D_{2})}\leq
C\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}$, it readily follows that

$R_{3}\leq
C\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}+C(\epsilon)\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\left\|\frac{\partial\theta_{v}(t)}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.$ (70)

To bound $R_{4}$ we notice that
$J_{h}-I_{h}^{(0)}J_{h}=(J_{h}-J)+(J-I_{h}^{(0)}J)+I_{h}^{(0)}(J-J_{h})$, so
by the triangle inequality it follows that

$\begin{array}[]{l}C\left\|(J_{h}-I_{h}^{(0)}J_{h})\theta_{vs}(t)\right\|_{L^{1}(D_{2})}\leq
C\left\|(J_{h}-J)\theta_{vs}(t)\right\|_{L^{1}(D_{2})}\\\ \\\
+C\left\|(J-I_{h}^{(0)}J)\theta_{vs}(t)\right\|_{L^{1}(D_{2})}+C\left\|I_{h}^{(0)}(J-J_{h})\theta_{vs}(t)\right\|_{L^{1}(D_{2})}.\end{array}$

Noting that
$\left\|ab\right\|_{L^{1}(D_{2})}\leq\frac{\varepsilon}{2}\left\|a\right\|_{L^{2}(D_{2})}^{2}+\frac{1}{2\varepsilon}\left\|b\right\|_{L^{2}(D_{2})}^{2}$
and applying the same argument as we have just done to bound $R_{3}$, we
obtain that

$\begin{array}[]{r}\left\|(J_{h}-I_{h}^{0}J_{h})\theta_{vs}(t)\right\|_{L^{1}(D_{2})}\leq
C\left(\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}+\left\|J-I_{h}^{0}J\right\|_{L^{2}(D_{2})}^{2}\right)\\\
\\\
+\sum_{l=1}^{M_{2}}\widehat{h}_{l}\theta_{v}^{2}(x_{l};R_{s}(x_{l}),t).\end{array}$

We bound the last term of this inequality as we have done for $R_{3},$ and by
virtue of assumption A3 set
$\left\|J-I_{h}^{(0)}J\right\|_{L^{2}(D_{2})}^{2}\leq
Ch^{2}\left\|\frac{\partial J}{\partial x}\right\|_{L^{2}(D_{2})}^{2}$. Hence,

$\begin{array}[]{r}R_{4}\leq
C(\epsilon)\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\displaystyle\left|\frac{\partial\theta_{v}(t)}{\partial
r}\right|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\\\ \\\
+C\left\|J-J_{h}\right\|_{L^{2}(D_{2})}^{2}+Ch^{2}\displaystyle\left\|\frac{\partial
J}{\partial x}\right\|_{L^{2}(D_{2})}^{2}.\end{array}$ (71)

Letting $\epsilon=\underline{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 $\theta_{u}$. Thus, subtracting
(15) from (10) it readily follows that

$\int_{D_{1}}\frac{\partial\theta_{u}}{\partial{t}}w_{h}dx+\int_{D_{1}}k_{1}\frac{\partial\theta_{u}}{\partial
x}\frac{dw_{h}}{dx}dx=\lambda\int_{D_{1}}\rho_{u}w_{h}dx-\int_{D_{1}}\frac{\partial\rho_{u}}{\partial{t}}w_{h}dx+\int_{D_{1}}a_{1}\left(J-J_{h}\right)w_{h}dx.$

Setting $w_{h}=\theta_{u}$ in this equation yields

$\begin{array}[]{r}\displaystyle\frac{1}{2}\frac{d}{dt}\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\underline{k_{1}}\displaystyle\left\|\frac{\partial\theta_{u}(t)}{\partial
x}\right\|_{L^{2}(D_{1})}^{2}\leq
C\displaystyle\left(\left\|\rho_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|\frac{\partial\rho_{u}(t)}{\partial
t}\right\|_{L^{2}(D_{1})}^{2}\right)\\\ \\\
+C\left\|J-J_{h}\right\|_{L^{2}(D_{1})}^{2}+C\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}.\end{array}$

By virtue of Lemma 11 we have the following result.

###### Lemma 13

There exists a constant $C$ independent of $h$ and $\Delta r$, but depending
on $\underline{k_{1}}$ and $\underline{k_{2}}$, such that

$\begin{array}[]{l}\displaystyle\frac{d}{dt}\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\displaystyle\underline{k_{1}}\left\|\frac{\partial\theta_{u}(t)}{\partial
x}\right\|_{L^{2}(D_{1})}^{2}\\\ \\\ \leq
C\left\\{h^{2}(\left\|\phi_{1}(t)\right\|_{H^{2}(D_{1})}^{2}+\left\|\phi_{2}(t)\right\|_{H^{2}(D_{2})}^{2}+\left|v_{s}(t)\right|_{H^{1}(D_{2})}^{2})\right.\\\
\\\
\left.+h^{4}\left(\left\|u(t)\right\|_{H^{2}(D_{1})}^{2}+\displaystyle\left\|\frac{\partial
u(t)}{\partial t}\right\|_{H^{2}(D_{1})}^{2}\right)+\Delta
r^{2}\left\|v(t)\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right\\}\\\
\\\
+C\left(\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\right)+\displaystyle\frac{\underline{k_{2}}}{4}\left\|\frac{\partial{\theta_{v}(t)}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}.\end{array}$ (72)

We are now in a position to establish the main result of this subsection.

###### Theorem 14

Let $(u,v,\phi_{1},\phi_{2})$ and $(u_{h},v_{h\Delta r},\phi_{1h},\phi_{2h})$
be the solutions to (10)-(13) and (15)-(18) respectively, with

$\left\|u(0)-u_{h}(0)\right\|_{L^{2}(D_{1})}\leq Ch^{2}\ \mathrm{and\
}\left\|v(0)-v_{h\Delta
r}(0)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\leq C(h+\Delta r^{2}).$

Furthermore, for $0\leq t\leq T_{\mathrm{end}},$ the following regularity
assumptions hold:

R1)$\ u$ and $\displaystyle\frac{\partial u}{\partial t}\in
L^{2}(0,T_{\mathrm{end}};H^{2}(D_{1})),$

R2) $v$ and $\displaystyle\frac{\partial v(t)}{\partial t}\in
L^{2}(0,T_{\mathrm{end}};L^{2}(D_{2};H_{r}^{2}(0,R_{s}(\cdot)))),$ and
$v_{s}\in L^{2}(0,T_{\mathrm{end}};H^{1}(D_{2})),$

R3) $\phi_{1}\in L^{2}(0,T_{\mathrm{end}};H^{2}(D_{1}))$, $\phi_{2}\in
L^{2}(0,T_{\mathrm{end}};H^{2}(D_{2})),$ and $J\in
L^{2}(0,T_{\mathrm{end}};H^{1}(D_{2}));$

then there is a constant $C(t,\underline{k_{1}},\underline{k_{2}})$ such that

$\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|v(t)-v_{h\Delta
r}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\int_{0}^{t}\left\|\Phi(\tau)-\Phi_{h}(\tau)\right\|_{V}^{2}d\tau\leq
C(h^{2}+\Delta r^{2}).$ (73)

Proof. Bounding $\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}$ as
$2\left(\left\|\rho_{u}(t)\right\|^{2}+\left\|\theta_{u}(t)\right\|^{2}\right)$
and using the estimate of Lemma 11 for
$\left\|v_{s}(t)-v_{sh}(t)\right\|_{L^{2}(D_{2})}$, Theorem 10 yields

$\begin{array}[]{l}\left\|\Phi(t)-\Phi_{h}(t)\right\|_{V}^{2}\leq
Ch^{2}\left|v_{s}(t)\right|_{H^{1}(D_{2})}^{2}+Ch^{4}\left\|u(t)\right\|_{H^{2}(D_{1})}^{2}+C\left\|\theta_{u}\right\|_{L^{2}(D_{1})}^{2}\\\
\\\
+C(\epsilon)\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\displaystyle\left\|\frac{\partial\theta_{v}(t)}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.\end{array}$

Considering this estimate together with those of Lemmas 12 and 13 it readily
follows that

$\begin{array}[]{r}\displaystyle\frac{d}{dt}\left(\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right)+\left\|\Phi(t)-\Phi_{h}(t)\right\|_{V}^{2}\leq
C(h^{2}+\Delta r^{2})\\\ \\\
+C\left(\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right).\end{array}$

Then, Gronwall inequality yields

$\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\int_{0}^{t}\left\|\Phi(\tau)-\Phi_{h}(\tau)\right\|_{V}^{2}d\tau\leq
C(t,\underline{k_{1}},\underline{k_{2}})(h^{2}+\Delta r^{2}).$ (74)

The terms $\left\|\theta_{u}(0)\right\|_{L^{2}(D_{1})}^{2}$ and
$\left\|\theta_{v}(0)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}$
are considered to be zero. From (57), (58), (59) and (35) it follows that

$\begin{array}[]{r}\left\|u(t)-u_{h}(t)\right\|_{L^{2}(D_{1})}^{2}+\left\|v(t)-v_{h\Delta
r}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\leq
2\left(\left\|\rho_{u}(t)\right\|_{L^{2}(D_{1})}+\left\|\theta_{u}(t)\right\|_{L^{2}(D_{1})}\right)\\\
\\\
+C\left(h^{2}\left\|v(t)\right\|_{H^{1}(D_{2};L_{r}^{2}(0,R_{s}(\cdot)))}+\left\|\rho_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot))}^{2}+\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot))}^{2}\right).\end{array}$

So, by combining this inequality with (74) we obtain (73).

##  5 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\Delta t$,
$\Delta 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,\ldots,N$, the
solution $(u_{h}^{n-1},v_{h\Delta r}^{n-1},\phi_{1h}^{n-1},\phi_{1h}^{n-1})\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})\times
W_{h}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$ is known,
calculate $(u_{h}^{n},v_{h\Delta r}^{n},\phi_{1h}^{n},\phi_{1h}^{n})\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})\times
W_{h}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$ as solution of
the system

$\int_{D_{1}}\widetilde{\partial}_{t}u^{n}{{}_{h}}w_{h}dx+\int_{D_{1}}k_{1}\frac{du_{h}^{n}}{dx}\frac{dw_{h}}{dx}dx=\int_{D_{1}}a_{1}J_{h}^{n}w_{h}dx\
\forall w_{h}\in V_{h}^{(1)}(\overline{D}_{1}).$ (75)

$\left\\{\begin{array}[]{l}\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\displaystyle\widetilde{\partial}_{t}v^{n}_{h\Delta
r}w_{h\Delta
r}r^{2}drdx+\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}k_{2}\displaystyle\frac{\partial
v^{n}{{}_{h\Delta r}}}{\partial{r}}\frac{\partial w_{h\Delta r}}{\partial
r}r^{2}drdx\\\ \\\
=-\displaystyle\int_{D_{2}}\frac{R_{s}^{2}(x)J_{h}^{n}w_{h\Delta
r}(x,R_{s}(x))}{a_{2}(x)F}dx\ \ \forall w_{h\Delta r}\in V_{h\Delta
r}(\overline{D}_{3}).\end{array}\right.$ (76)

$\int_{D_{1}}\kappa(u_{h}^{n})\frac{d\phi_{1h}^{n}}{dx}\frac{dw_{h}}{dx}dx=\int_{D_{1}}J_{h}^{n}w_{h}dx\
\ \forall w_{h}\in V_{h}^{(1)}(\overline{D}_{1}).$ (77)

$\int_{D_{2}}\sigma\frac{d\phi_{2h}^{n}}{dx}\frac{dw_{h}}{dx}dx=-\int_{D_{2}}\left(J_{h}^{n}+g\right)w_{h}dx\
\ \forall w_{h}\in V_{h}^{(1)}(\overline{D}_{2}).$ (78)

$\int_{D_{2}}J_{h}^{n}dx=0\text{ with
}\int_{D_{\mathrm{a}}}J_{h}^{n}dx=I(t_{n})=-\int_{D_{\text{{c}}}}J_{h}^{n}dx,$
(79)

where

$\begin{array}[]{l}\widetilde{\partial}_{t}u_{h}^{n}=\displaystyle\frac{u_{h}^{n}-u_{h}^{n-1}}{\Delta
t},\text{{\ \ }}\widetilde{\partial}_{t}v_{h\Delta
r}^{n}=\displaystyle\frac{v_{h\Delta r}^{n}-v_{h\Delta r}^{n-1}}{\Delta t},\\\
\\\
J_{h}^{n}=J(x,u_{{}_{h}}^{n},v_{sh}^{n},\eta_{h}^{n})=a_{2}(x)i_{0h}^{n}\sinh\left(\beta\eta_{h}^{n}\right),\
i_{0h}^{n}=i_{0}(u_{h}^{n},v_{sh}^{n}),\\\ \\\
\eta_{h}^{n}=\phi_{1h}^{n}-\phi_{2h}^{n}-\alpha^{n}_{h}\ln
u_{h}^{n}-\overline{U}_{h}(v_{sh}^{n}),\
\alpha^{n}_{h}=\alpha(u_{h}^{n}).\end{array}$ (80)

###  5.1 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_{h}^{n},v_{h\Delta r}^{n},v_{sh}^{n})\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})\times
V_{h}^{(0)}(\overline{D}_{2})$ and the assumptions A1-A4 hold, the system
(77)-(78) has a unique solution $(\phi_{1h}^{n},\phi_{2h}^{n})\in
W_{h}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$; then, returning
to the system (75)-(76) and applying a well-known consequence of Brower´s
fixed point theorem, which is presented as Corollary 1.1 in [8], we prove that
there exists $(u_{h}^{n},v_{h\Delta r}^{n})\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})$.

###### Lemma 15

Assuming that for all $n$, $(u_{h}^{n},v_{h\Delta r}^{n},v_{sh}^{n})\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})\times
V_{h}^{(0)}(\overline{D}_{2})$, and the assumptions A1-A4 hold, then the
system (77)-(78) has a unique solution $(\phi_{1h}^{n},\phi_{2h}^{n})\in
W_{h}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$.

Proof. Looking at (77)-(78) and in order to apply Minty-Browder theorem to
prove the existence of a solution, we define the functions

$\widehat{\eta}_{h}^{n}:=-\alpha_{h}^{n}\ln
u_{h}^{n}-\overline{U}_{h}(v_{sh}^{n})\ \mathrm{and\
}\widehat{J}_{h}^{n}:=J_{h}^{n}(x,u_{h}^{n},v_{sh}^{n},\widehat{\eta}_{h}^{n}),$

it is worth noticing that $\widehat{\eta}_{h}^{n}$ is equal to $\eta_{h}^{n}$
when the potentials $\phi_{1h}^{n}$ and $\phi_{1h}^{n}$ are zero. Now, going
back to Section 4.2 and using $\widehat{J}_{h}^{n}$, we define the operators
$\widehat{B}_{h}:V_{h}\rightarrow V_{h}^{\ast}$ and $A_{h}:V_{h}\rightarrow
V_{h}^{\ast}$ as follows: for all $n=1,2,..,N$,

$\left\langle\widehat{B}_{h}(\Phi_{h}^{n}),\Psi_{h}\right\rangle=\int_{D_{2}}(J_{h}^{n}-\widehat{J}_{h}^{n})(\psi_{2h}-\psi_{1h})dx\
\ \forall\Psi_{h}\in V_{h}.$

and

$\left\langle
A_{h}(\Phi_{h}^{n}),\Psi_{h}\right\rangle=a_{h}(\Phi_{h}^{n},\Psi_{h})+\left\langle\widehat{B}_{h}(\Phi_{h}^{n}),\Psi_{h}\right\rangle.$

Notice that when $\Phi_{h}^{n}=(0,0)$,
$\left\langle\widehat{B}_{h}(\Phi_{h}^{n}),\Psi_{h}\right\rangle=0$ because
$J_{h}^{n}=\widehat{J}_{h}^{n}$. Now, we can recast (77)-(78) as follows. Find
$\Phi_{h}^{n}:=(\phi_{1}^{n},\phi_{1}^{n})\in W_{h}(\overline{D}_{1})\times
V_{h}^{(1)}(\overline{D}_{2})$ such that

$\left\langle
A_{h}(\Phi_{h}^{n}),\Psi_{h}\right\rangle=-\int_{D_{2}}g\psi_{2h}dx-\int_{D_{2}}\widehat{J}_{h}^{n}(\psi_{2h}-\psi_{1h})dx\
\ \forall\Psi_{h}\in V_{h}.$ (81)

We can prove, using the same arguments as in Lemma 8, that the operator
$\widehat{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.,
$\forall\Phi_{h}^{n}\in W_{h}(\overline{D}_{1})\times
V_{h}^{(1)}(\overline{D}_{2})$, there exists a positive constant $\alpha$ such
that

$\left\langle
A_{h}(\Phi_{h}^{n}),\Phi_{h}^{n}\right\rangle\geq\alpha\left\|\Phi_{h}^{n}\right\|_{V}^{2}.$

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 $\forall\Phi_{h}^{n}$, $\overline{\Phi}_{h}^{n}$ $\in
W_{h}(\overline{D}_{1})\times V_{h}^{(1)}(\overline{D}_{2})$

$a_{h}(\Phi_{h}^{n}-\overline{\Phi}_{h}^{n},\Phi_{h}^{n}-\overline{\Phi}_{h}^{n})+\left\langle\widehat{B}_{h}(\Phi_{h}^{n})-\widehat{B}_{h}(\overline{\Phi}_{h}^{n}),\Phi_{h}^{n}-\overline{\Phi}_{h}^{n}\right\rangle\geq\alpha\left\|\Phi_{h}^{n}-\overline{\Phi}_{h}^{n}\right\|_{V}^{2},$

then taking$\ \overline{\Phi}_{h}^{n}=(0,0)$ it follows the coerciveness of
$A_{h}$. Hence, the Minty-Browder theorem [18] guaranties the existence of a
solution $\Phi_{h}^{n}$ 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
$\Phi_{h}^{n}:=(\phi_{1h}^{n},\phi_{2h}^{n})$ and
$\overline{\Phi}_{h}^{n}:=\left(\overline{\phi}_{1h}^{n},\overline{\phi}_{2h}^{n}\right)$
of (77)-(78), then setting, $z_{1h}=\phi_{1h}^{n}-\overline{\phi}_{1h}^{n}$
and $z_{2h}=\phi_{2h}^{n}-\overline{\phi}_{2h}^{n}$, from (77) it follows that

$\int_{D_{1}}\kappa(u_{h}^{n})\frac{dz_{1h}}{dx}\frac{dw_{h}}{dx}dx=\int_{D_{2}}\left(J_{h}^{n}-\overline{J}_{h}^{n}\right)w_{h}dx\
\ \forall w_{h}\in V_{h}^{\left(1\right)}(\overline{D}_{1})$

and from (78)

$\int_{D_{1}}\sigma\frac{dz_{2h}}{dx}\frac{dv_{h}}{dx}dx=-\int_{D_{2}}\left(J_{h}^{n}-\overline{J}_{h}^{n}\right)v_{h}dx\
\ \forall v_{h}\in V_{h}^{\left(1\right)}(\overline{D}_{2}),$

where $J_{h}^{n}=J_{h}^{n}(x,u_{h}^{n},v_{sh}^{n},\eta_{h}^{n})$ and
$\overline{J}_{h}^{n}=J_{h}^{n}(x,u_{h}^{n},v_{sh}^{n},\overline{\eta}_{h}^{n})$,
with
$\overline{\eta}_{h}^{n}=\overline{\phi}_{2h}^{n}-\overline{\phi}_{1h}^{n}-\alpha_{h}^{n}\ln
u_{h}^{n}-\overline{U}_{h}(v_{sh}^{n})$. Setting $w_{h}=z_{1h}$ and
$v_{h}=z_{2h}$ and applying the mean value theorem one readily obtains that

$\int_{D_{1}}\kappa(u_{h}^{n})\left(\frac{dz_{1h}}{dx}\right)^{2}dx+\int_{D_{1}}\sigma\left(\frac{dz_{2h}}{dx}\right)^{2}dx+\int_{D_{2}}\frac{\partial
J_{h}^{n}\left(\xi\right)}{\partial\eta^{n}}(z_{2h}-z_{1h})^{2}=0,$

here $\frac{\partial J_{h}^{n}\left(\xi\right)}{\partial\eta^{n}}>0$ according
to assumption A3. The first term of this expression implies that for all $n$,

$z_{1h}=\phi_{1h}^{n}-\overline{\phi}_{1h}^{n}=K_{1},$

but the constant $K_{1}=0$ because $\phi_{1h}^{n}$ and
$\overline{\phi}_{1h}^{n}$ are in $W_{h}(\overline{D}_{1})$, so
$\phi_{1h}^{n}=\overline{\phi}_{1h}^{n}$ Similarly, from the second and third
terms it follows that $z_{2h}=0$, and consequently
$\phi_{2h}^{n}=\overline{\phi}_{2h}^{n}$. Hence, we have just proved that for
all $n$ there is a unique solution $(\phi_{1h}^{n},\phi_{2h}^{n})$.

###### Lemma 16

Let $(\phi_{1h}^{n},\phi_{2h}^{n})\in W_{h}(\overline{D}_{1})\times
V_{h}^{(1)}(\overline{D}_{2})$ be the solution to (77)-(78). There exists a
unique solution $\left(u_{h}^{n},v_{h\Delta r}^{n}\right)\in
V_{h}^{(1)}(\overline{D}_{1})\times V_{h\Delta r}(\overline{D}_{3})$ to the
system (75)-(76).

Proof. We start proving the existence of $u_{h}^{n}\in
V_{h}^{(1)}(\overline{D}_{1})$ as solution of (75). To this end, we write (75)
as $F_{h}(u_{h}^{n})=0$, where $F_{h}:V_{h}^{(1)}(\overline{D}_{1})\rightarrow
V_{h}^{(1)}(\overline{D}_{1})$ is a continuous mapping defined by the relation

$\begin{array}[]{r}\displaystyle\int_{D_{1}}F_{h}(\chi_{h})w_{h}dx=\displaystyle\int_{D_{1}}\left(\chi_{h}-u_{h}^{n-1}\right)w_{h}dx+\Delta
t\int_{D_{1}}k_{1}\frac{d\chi_{h}^{n}}{dx}\frac{dw_{h}}{dx}dx\\\ \\\
-\displaystyle\Delta t\int_{D_{1}}a_{1}J_{h}^{n}(\chi_{h})w_{h}dx=0\ \ \forall
w_{h}\in V_{h}^{(1)}(\overline{D}_{1}),\end{array}$

here,
$J_{h}^{n}(\chi_{h})=J(x,\chi_{h},v_{sh},\phi_{1h}^{n},\phi_{2h}^{n},\overline{U}_{h}(v_{sh})),$
with $v_{sh}$ being picked up from $S_{Q}^{\ast}$ because we assume that
$v_{h\Delta r}^{n}$ belongs to this space; moreover, we also assume that
$\chi_{h}$ is in $S_{P}$. According to Brower´s fixed point theorem, the
equation $F_{h}(\chi_{h})=0$ has a solution$\ \chi_{h}\in
B_{q}:=\left\\{v_{h}\in
V_{h}^{(1)}(\overline{D}_{1}):\left\|v_{h}\right\|_{L^{2}(D_{1})}\leq
q\right\\}$, if $\int_{D_{1}}F_{h}(\chi_{h})\chi_{h}dx>0$ for
$\left\|\chi_{h}\right\|_{L^{2}(D_{1})}=q$. On account of the assumptions
$v_{sh}\in S_{Q}^{\ast}$ and $\chi_{h}\in S_{P}$, it follows that there exists
a constant $C_{1}=C_{1}(P,Q,K)$ such that
$\int_{D_{1}}a_{1}J_{h}^{n}(\chi_{h})\chi_{h}dx\leq
C_{1}\left\|\chi_{h}\right\|_{L^{2}(D_{1})}$. Hence,

$\begin{array}[]{r}\displaystyle\int_{D_{1}}F_{h}(\chi_{h})\chi_{h}dx\geq\left\|\chi_{h}\right\|_{L^{2}(D_{1})}^{2}-\left\|u_{h}^{n-1}\right\|_{L^{2}(D_{1})}^{2}+\Delta
tk_{0}\left\|\frac{d\chi_{h}^{n}}{dx}\right\|_{L^{2}(D_{1})}^{2}-\Delta
tC_{1}\left\|\chi_{h}\right\|_{L^{2}(D_{1})}\\\ \\\
\geq\left\|\chi_{h}\right\|_{L^{2}(D_{1})}^{2}-\left\|u_{h}^{n-1}\right\|_{L^{2}(D_{1})}^{2}-\Delta
tC_{1}\left(1+\left\|\chi_{h}\right\|_{L^{2}(D_{1})}\right)\left\|\chi_{h}\right\|_{L^{2}(D_{1})}.\end{array}$

Then, taking $\Delta t\leq\Delta t_{0}<1/C_{1}$,
$\int_{D_{1}}F_{h}(\chi_{h})\chi_{h}dx$ is positive for
$\left\|\chi_{h}\right\|_{L^{2}(D_{1})}$ sufficiently large. This shows the
existence of the solution $u_{h}^{n}\in V_{h}^{(1)}(\overline{D}_{1})$. Next,
we prove the uniqueness. To this end, we consider that there exist $X$ and $Y$
$\in V_{h}^{(1)}(\overline{D}_{1})$ satisfying (75), so

$\int_{D_{1}}\left(X-Y\right)w_{h}dx+\Delta
t\int_{D_{1}}k_{1}\frac{d(X-Y)}{dx}\frac{dw_{h}}{dx}dx=\Delta
t\int_{D_{1}}a_{1}\left(J_{h}^{n}(X)-J_{h}^{n}(Y)\right)w_{h}dx\ \ \forall
w_{h}\in V_{h}^{(1)}(\overline{D}_{1}).$

Setting $w_{h}=X-Y$ and invoking the arguments of Lemmas 5 and 6 yields

$\left\|X-Y\right\|_{L^{2}(D_{1})}^{2}+k_{1}\Delta
t\left\|\frac{d(X-Y)}{dx}\right\|_{L^{2}(D_{1})}^{2}\leq\Delta
tC_{2}\left\|X-Y\right\|_{L^{2}(D_{1})}^{2},$

where the constant $C_{2}=C_{2}(P,Q,K)$. Thus, taking $\Delta t\leq\Delta
t_{0}<1/C_{2}$ it follows that $X=Y$. It remains to prove the existence and
uniqueness of $v_{h\Delta r}^{n}$, but the arguments to be used for such a
proof are the same as for $u_{h}^{n}$, so we omit them.

###  5.2 Error estimates for the fully discrete solution

As in Section 4.3, we write for $t=t_{n}$

$\left\\{\begin{array}[]{l}u^{n}-u_{h}^{n}=\rho_{u}^{n}+\theta_{u}^{n},\\\ \\\
v^{n}-v_{h\Delta
r}^{n}=v^{n}-I_{0}^{x}v^{n}+I_{0}^{x}\rho_{v}^{n}+\theta_{v}^{n}.\end{array}\right.$
(82)

###### Theorem 17

Let $(u_{h}^{n},v_{h\Delta r}^{n},\phi_{1}^{n},\phi_{2}^{n})$ be the solution
to (75)-(78). Then, under proper regularity assumptions there exists a
constant $C$ such that for $\Delta t$ small

$\left\|u^{n}-u_{h}^{n}\right\|_{L^{2}(D_{1})}^{2}+\left\|v^{n}-v_{h\Delta
r}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\Delta
t\sum_{j=1}^{t_{n}}\left\|\Phi^{j}-\Phi_{h}^{j}\right\|_{V}^{2}d\tau\leq
C(h^{2}+\Delta r^{2}+\Delta t^{2}).$ (83)

The constant $C$ is of the form
$C(\Gamma)\exp(C(\underline{k_{1}},\underline{k_{2}})t_{n}$, $C(\Gamma)$ being
another constant that depends on the exact solution $(u,v,\phi_{1},\phi_{2})$,
see (87) below.

Proof. Since $\rho_{u}^{n}$, $v^{n}-I_{0}^{x}v^{n}+I_{0}^{x}\rho_{v}^{n}$ and
$\left\|\Phi^{j}-\Phi_{h}^{j}\right\|_{V}^{2}$ are estimated as in Section
4.3, we shall address our attention to the estimations for $\theta_{u}^{n}$
and $\theta_{v}^{n}$. We start with the calculation for $\theta_{v}^{n}$. For
this purpose, we recast (64) for $t=t_{n}$ as

$\begin{array}[]{c}\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(\widetilde{\partial}_{t}I_{0}^{x}{v}^{n}w+k_{2}\frac{\partial
I_{0}^{x}v^{n}}{\partial r}\frac{\partial w}{\partial r}\right)r^{2}drdx\\\
\\\
=-\displaystyle\int_{D_{2}}I_{0}^{x}(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}J^{n}w(x;R_{s}(x)))dx\\\
\\\
+\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(\widetilde{\partial}_{t}I_{0}^{x}{v}^{n}w-\frac{\partial
I_{0}^{x}v^{n}}{\partial t}w\right)r^{2}drdx.\end{array}$ (84)

Setting, as we did in Section 4.3, $v_{h\Delta
r}^{n}=I_{0}^{x}v^{n}-\left(I_{0}^{x}\rho_{v}^{n}+\theta_{v}^{n}\right)$ in
(76) and using (32) and (84 ) yields for $t=t_{n}$

$\begin{array}[]{r}\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(\widetilde{\partial}_{t}\theta_{v}^{n}w_{h\Delta
r}+k_{2}\frac{\partial\theta_{v}^{n}}{\partial r}\frac{\partial w_{h\Delta
r}}{\partial
r}\right)r^{2}drdx=\lambda\int_{D_{2}}\int_{0}^{R_{s}(x)}I_{0}^{x}\rho_{v}^{n}w_{h\Delta
r}r^{2}drdx\\\ \\\
-\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\widetilde{\partial}_{t}I_{0}^{x}\rho_{v}^{n}w_{h\Delta
r}r^{2}drdx\\\ \\\
-\displaystyle\int_{D_{2}}I_{h}^{0}(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}(J^{n}(x)-J_{h}^{n}(x))w_{hs}(x))dx\\\
\\\
+\displaystyle\int_{D_{2}}\left(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}J_{h}^{n}(x)-I_{h}^{0}(R_{s}^{2}(x)a_{2}^{-1}(x)F^{-1}J_{h}^{n}(x))\right)w_{hs}(x)dx,\\\
\\\
-\displaystyle\int_{D_{2}}\int_{0}^{R_{s}(x)}\left(\widetilde{\partial}_{t}I_{0}^{x}v^{n}-\frac{\partial
I_{0}^{x}v^{n}}{\partial t}\right)r^{2}drdx\end{array}$

Letting $w_{h\Delta r}=\theta_{v}^{n}$, $w_{hs}=\theta_{vs}^{n}$, and noting
that for $a$ and $b$ real numbers,
$2(a-b)b=a^{2}-b^{2}{\color[rgb]{0,0,1}-}(a-b)^{2}$, it follows that

$\begin{array}[]{l}\displaystyle\frac{1}{2}\widetilde{\partial}_{t}\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\displaystyle\underline{k_{2}}\left\|\frac{\partial\theta_{v}^{n}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\\\ \\\
\leq\lambda\left\|I_{0}^{x}\rho_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\\\
\\\
+\displaystyle\left\|\widetilde{\partial}_{t}I_{0}^{x}\rho_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\\\
\\\
+C\left\|I_{h}^{0}\left((J^{n}-J_{h}^{n})\theta_{vs}^{n}\right)\right\|_{L^{1}(D_{2})}+C\left\|(J_{h}^{n}-I_{h}^{0}J_{h}^{n})\theta_{vs}^{n}\right\|_{L^{1}(D_{2})}\\\
\\\ +\left\|\widetilde{\partial}_{t}I_{0}^{x}v^{n}-\frac{\partial
I_{0}^{x}v^{n}}{\partial{t}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\equiv\sum_{i=1}^{5}R_{i}^{n}.\end{array}$
(85)

We bound the right hand side of this inequality applying the same arguments as
in (66). Thus, we have that

$R_{1}^{n}\leq C\Delta
r^{4}\left\|v^{n}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}+C\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}$

$R_{2}^{n}\leq\frac{C\Delta r^{4}}{\Delta
t}\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial v}{\partial
t}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}dt+C\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.$

$R_{3}^{n}\leq
C\left\|J^{n}-J_{h}^{n}\right\|_{L^{2}(D_{2})}^{2}+C(\epsilon)\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\left\|\frac{\partial\theta_{v}^{n}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.$

$\begin{array}[]{r}R_{4}^{n}\leq
C\left\|J^{n}-J_{h}^{n}\right\|_{L^{2}(D_{2})}^{2}+C(\epsilon)\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\displaystyle\left|\frac{\partial\theta_{v}^{n}}{\partial
r}\right|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\\\ \\\
+Ch^{2}\displaystyle\left\|\frac{\partial J^{n}}{\partial
x}\right\|_{L^{2}(D_{2})}^{2}.\end{array}$

In both $R_{3}^{n}$ and $R_{4}^{n}$ the term
$C\left\|J^{n}-J_{h}^{n}\right\|_{L^{2}(D_{2})}^{2}$ is bounded by Lemma 11
for $t=t_{n}$; thus, using the notation
$\Gamma=(u,v,v_{s},\phi_{1},\phi_{2})$, we can set that

$\begin{array}[]{l}C\left\|J^{n}-J_{h}^{n}\right\|_{L^{2}(D_{2})}^{2}\leq
C(\Gamma)(h^{2}+\Delta
r^{2})+C\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}\\\ \\\
+C(\epsilon)\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\displaystyle\left\|\frac{\partial\theta_{v}^{n}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2},\end{array}$ (86)

where the constant $C(\Gamma)$ is given as

$\begin{array}[]{c}C(\Gamma)=C\max\left(\left\|\phi_{1}\right\|_{L^{\infty}(0,T_{\mathrm{end}};H^{2}(D_{1}))},\left\|\phi_{2}\right\|_{L^{\infty}(0,T_{\mathrm{end}};H^{2}(D_{2}))},\left\|v_{s}\right\|_{L^{\infty}(0,T_{\mathrm{end}};H^{1}(D_{2}))},\right.\\\
\\\
\left.\left\|u\right\|_{L^{\infty}(0,T_{\mathrm{end}};H^{2}(D_{1}))},\left\|v\right\|_{L^{\infty}(0,T_{\mathrm{end}};L^{2}(D_{2};H_{r}^{2}(0,R_{s}(\cdot))))}\right).\end{array}$
(87)

Hence, we can write

$\begin{array}[]{r}R_{3}^{n}+R_{4}^{n}\leq C(\Gamma)(h^{2}+\Delta
r^{2})+Ch^{2}\displaystyle\left\|\frac{\partial J^{n}}{\partial
x}\right\|_{L^{2}(D_{2})}^{2}+C\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}\\\
\\\
+C(\epsilon)\left\|\theta_{v}(t)\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\epsilon\displaystyle\left\|\frac{\partial\theta_{v}(t)}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.\end{array}$

To estimate the term $R_{5}^{n}$, we notice that by approximation theory

$\left\|\widetilde{\partial}_{t}I_{0}^{x}v^{n}-\frac{\partial
I_{0}^{x}v^{n}}{\partial{t}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\leq
C\left\|\widetilde{\partial}_{t}v^{n}-\frac{\partial
v^{n}}{\partial{t}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}$

and

$\widetilde{\partial}_{t}v^{n}-\frac{\partial
v^{n}}{\partial{t}}=\frac{-1}{\Delta
t}\int_{t_{n-1}}^{t_{n}}(t-t_{n-1})\frac{\partial^{2}v}{\partial t^{2}}dt,$

so,

$\left\|\widetilde{\partial}_{t}I_{0}^{x}v^{n}-\frac{\partial
I_{0}^{x}v^{n}}{\partial{t}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}\leq
C\left(\Delta
t\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial^{2}v}{\partial{t}^{2}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}dt\right)^{1/2}.$

Applying Young inequality yields

$R_{5}^{n}\leq C\Delta
t\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial^{2}v}{\partial{t}^{2}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}dt+C\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}$

Collecting these bounds in (85) and letting $\epsilon=\underline{k_{2}}/2$
yields

$\begin{array}[]{r}\displaystyle\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\displaystyle\Delta
t\underline{k_{2}}\left\|\frac{\partial\theta_{v}^{n}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\leq\displaystyle\left\|\theta_{v}^{n-1}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+F_{v}^{n}\\\
\\\ +C(\underline{k_{2}})\Delta
t\left(\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}\right),\end{array}$
(88)

where

$\begin{array}[]{c}F_{v}^{n}=\Delta t\left(C(\Gamma)(h^{2}+\Delta
r^{2})+Ch^{2}\displaystyle\left\|\frac{\partial J^{n}}{\partial
x}\right\|_{L^{2}(D_{2})}^{2}+C\Delta
r^{4}\left\|v^{n}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right)\\\
\\\ +C\Delta r^{4}\displaystyle\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial
v}{\partial t}\right\|_{L^{2}(D_{2},H_{r}^{2}(0,R_{s}(\cdot)))}^{2}dt+C\Delta
t^{2}\displaystyle\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial^{2}v}{\partial{t}^{2}}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}dt.\end{array}$
(89)

To calculate an estimate for $\theta_{u}^{n}$, we observe that subtracting
(75) from (10) and setting $\theta_{u}^{n}=e_{u}^{n}-\rho_{u}^{n}$ it follows
that

$\begin{array}[]{r}\displaystyle\int_{D_{1}}\widetilde{\partial}_{t}\theta_{u}^{n}w_{h}dx+\displaystyle\int_{D_{1}}k_{1}\frac{d\theta_{u}^{n}}{dx}\frac{dw_{h}}{dx}dx=\lambda\int_{D_{1}}\rho_{u}^{n}w_{h}dx-\int_{D_{1}}\partial\rho_{u}^{n}w_{h}dx\\\
\\\
+\displaystyle\int_{D_{1}}a_{1}\left(J^{n}-J_{h}^{n}\right)w_{h}dx+\displaystyle\int_{D_{1}}\left(\widetilde{\partial}_{t}u^{n}-\frac{\partial
u^{n}}{\partial t}\right)w_{h}dx.\end{array}$

Letting $w_{h}=\theta_{u}^{n}$ yields

$\begin{array}[]{r}\displaystyle\frac{1}{2}\widetilde{\partial}_{t}\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}+\underline{k_{1}}\displaystyle\left\|\frac{d\theta_{u}^{n}}{dx}\right\|_{L^{2}(D_{1})}^{2}\leq
C\displaystyle\left(\left\|\rho_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}+\frac{1}{\Delta
t}\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial\rho_{u}}{\partial
t}\right\|_{L^{2}(D_{1})}^{2}dt\right)\\\ \\\ +C\Delta
t\displaystyle\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial^{2}u}{\partial^{2}t}\right\|_{L^{2}(D_{1})}^{2}dt+C\left\|J^{n}-J_{h}^{n}\right\|_{L^{2}(D_{1})}^{2}+C\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}.\end{array}$

Then by virtue of (29) and Lemma 11 it follows that

$\begin{array}[]{l}\left\|\theta_{u}^{n}(t)\right\|_{L^{2}(D_{1})}^{2}+\Delta
t\underline{k_{1}}\displaystyle\left\|\frac{d\theta_{u}^{n}}{dx}\right\|_{L^{2}(D_{1})}^{2}\leq\left\|\theta_{u}^{n-1}(t)\right\|_{L^{2}(D_{1})}^{2}+F_{u}^{n}\\\
\\\ +C(\underline{k_{1}},\underline{k_{2}})\Delta
t\left(\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}\right)+\displaystyle\frac{\underline{k_{2}}}{2}\left\|\frac{\partial{\theta_{v}^{n}}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\end{array}$ (90)

where

$F_{u}^{n}=\Delta tC(\Gamma)(h^{2}+\Delta
r^{2})+Ch^{4}\displaystyle\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial
u}{\partial t}\right\|_{H^{2}(D_{1})}^{2}dt+C\Delta
t^{2}\displaystyle\int_{t_{n-1}}^{t_{n}}\left\|\frac{\partial^{2}u}{\partial^{2}t}\right\|_{L^{2}(D_{1})}^{2}dt$
(91)

It remains to estimate $\left\|\Phi^{n}-\Phi_{h}^{n}\right\|_{V}$. Returning
to the proof of Theorem 14 we have that for $t=t_{n}$

$\left\|\Phi^{n}-\Phi_{h}^{n}\right\|_{V}^{2}\leq
C(\Gamma)h^{2}+C(\epsilon)\left(\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right)+\epsilon\displaystyle\left\|\frac{\partial\theta_{v}^{n}}{\partial
r}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}.$ (92)

Thus, setting $\epsilon=\underline{k_{2}}/2$ in (92) and adding (88), (90) and
(92) we obtain that

$\begin{array}[]{l}\left(1-C(\underline{k_{1}},\underline{k_{2}})\Delta
t\right)\left(\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}+\displaystyle\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right)+\Delta
t\left\|\Phi^{n}-\Phi_{h}^{n}\right\|_{V}^{2}\leq F_{v}^{n}+F_{u}^{n}\\\ \\\
+\left\|\theta_{u}^{n-1}\right\|_{L^{2}(D_{1})}^{2}+\displaystyle\left\|\theta_{v}^{n-1}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\end{array}$

For $\Delta t$ small

$\begin{array}[]{r}\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}+\displaystyle\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\Delta
t\left\|\Phi^{n}-\Phi_{h}^{n}\right\|_{V}^{2}\leq
C\left(F_{v}^{n}+F_{u}^{n}\right)\\\ \\\
\left(1+C(\underline{k_{1}},\underline{k_{2}})\Delta
t\right)\left(\left\|\theta_{u}^{n-1}\right\|_{L^{2}(D_{1})}^{2}+\displaystyle\left\|\theta_{v}^{n-1}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}\right).\end{array}$

Hence, by repeated application and taking
$\left\|\theta_{u}^{0}\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}^{0}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}=0$,
it results that

$\left\|\theta_{u}^{n}\right\|_{L^{2}(D_{1})}^{2}+\left\|\theta_{v}^{n}\right\|_{L^{2}(D_{2},L_{r}^{2}(0,R_{s}(\cdot)))}^{2}+\Delta
t\sum_{j=1}^{n}\left\|\Phi^{j}-\Phi_{h}^{j}\right\|_{V}^{2}\leq
C\displaystyle\sum_{j=1}^{n}\frac{F_{v}^{j}+F_{u}^{j}}{\left(1+C(\underline{k_{1}},\underline{k_{2}}))\Delta
t\right)^{j-n}}.$

Noting that $\left(1+C(k_{1},k_{2})\Delta t\right)\leq e^{C(k_{1},k_{2})\Delta
t}$, then we can write

$\begin{array}[]{c}\displaystyle\sum_{j=1}^{n}\frac{F_{v}^{j}+F_{u}^{j}}{\left(1+C(\underline{k_{1}},\underline{k_{2}})\Delta
t\right)^{j-n}}\leq
e^{C(\underline{k_{1}},\underline{k_{2}})t_{n}}\sum_{j=1}^{n}F_{v}^{j}+F_{u}^{j}\
\ \left(\mathrm{by\ (\ref{eef6.1})\ and\ (\ref{eef7})}\right)\\\ \\\ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \leq
C(\Gamma)e^{C(\underline{k_{1}},\underline{k_{2}})t_{n}}\left(h^{2}+\Delta
r^{2}+\Delta t^{2}\right).\end{array}$

This completes the proof

##  6 Acknowledgements

This research has been partially funded by grant PGC-2018-097565-B100 of
Ministerio de Ciencia, Innovación y Universidades of Spain.

## References

  * [1] R. Bermejo, P. Galán del Sastre, An implicit-explicit Runge-Kutta-Chebyshev finite element method for the nonlinear Lithium-ion battery equations, Applied Mathematics and Computation 36 (2019) 398-420.
  * [2] P. G. Ciarlet, The Finite Element Method for Elliptic Prolems, North Holland, Amsterdand 1978.
  * [3] M. Doyle, T. F. Fuller, J. Newman, Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell, J. Electrochem. Soc. 140 (1993) 1526-1533.
  * [4] J.I.Díaz, D. Gómez-Castro, A.M. Ramos, On the well-posedness of a multiscale mathematical model for Lithium-ion batteries, Adv. Nonlinear Anal. 8 (2019) 1132–1157.
  * [5] K. Eriksson, V. Thomée, Galerkin methods for singular value problems in one space dimension, Math. Comp. 42 (1984) 345-367.
  * [6] Z. Farkas, I. Faragó, Á. Kriston, A. Pfrang, Improvement of the accuracy of multi-scale models of Li-ion batteries by applying operating splitting techniques, J. Computational and Applied Mathjematics 310 (2017) 59-79.
  * [7] D. Gilbert, N. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, Heidelberg 1977.
  * [8] V. Girault, P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer-Verlag, Heidelberg 1986.
  * [9] S. Kosch , J. Sturm , J. Schuster , G. Mulder , E. Ayerbe , A. Jossem , A computationally efficient multi-scale model for lithium-ion cells, J. Electrochem. Soc. 165 (2018) A2374–A2388.
  * [10] C. Kr$\ddot{\text{o}}$ner, A mathematical exploration of a PDE system for lithium-ion batteries, Ph.D. Thesis in Mathematics University of California, Berkeley (2016) 63 pp.
  * [11] S. Mazunder, J. Lu, Faster-than-real-time simulation of Lithium-ion batteries with full spatial and temperal resolution. International Journal of Electrochemistry, 268747 (2013) 10 pp. ID 268747.
  * [12] J. Newman, Dualfoil 5.1 Newman research group webpage. http://www.cchem.berkeley.edu/jsngrp/.
  * [13] J. Newman, K.E. Thomas-Alyea, Electrochemical Systems, 3rd ed. Wiley, Hoboken, NJ. 2004.
  * [14] P.W.C. Northrop, M. Pathak, D. Rife, S. De, S. Santhanagopalan, V.R. Subramanian, Efficient simulation and model reformulation of two-dimensional electrochemical thermal behavior of Lithium-ion batteries. J. Electrochem.. Soc., 162 (2015) A940-A951.
  * [15] G.L. Plett, Battery Management System vol.1, Artech House, Boston. 2015\.
  * [16] R. Schreiber, S. Eisenstat, Finite element methods for spherically symmetric elliptic equations. SIAM J. Numer. Anal., 18 (1981) 546-558.
  * [17] K. A. Smith, C-Y, Wang, Solid-state diffusion limitations on pulse operation of a Lithium-ion cell for hybrid electric vehicles. Journal of Power Sources, 161 (2006) 628-639.
  * [18] E. Zeidler, Nonlinear Functional Analysis and Its Applications IIB, Springer-Verlag, New York, 2013.
  * [19] J. Wu, J. Xu, H.Zou, On the well-posedness of a mathematical model for lithium-ion battery systems. Methods and Applications of Analysis, 13 (2006) 275-298.