# Error estimates for semi-discrete finite element approximations for a moving
boundary problem capturing the penetration of diffusants into rubber

Surendra Nepal 1,∗, Yosief Wondmagegne1, Adrian Muntean1
1 Department of Mathematics and Computer Science, Karlstad University, Sweden
<EMAIL_ADDRESS>
###### Abstract

We consider a moving boundary problem with kinetic condition that describes
the diffusion of solvent into rubber and study semi-discrete finite element
approximations of the corresponding weak solutions. We report on both a priori
and a posteriori error estimates for the mass concentration of the diffusants,
and respectively, for the a priori unknown position of the moving boundary.
Our working techniques include integral and energy-based estimates for a
nonlinear parabolic problem posed in a transformed fixed domain combined with
a suitable use of the interpolation-trace inequality to handle the interface
terms. Numerical illustrations of our FEM approximations are within the
experimental range and show good agreement with our theoretical investigation.
This work is a preliminary investigation necessary before extending the
current moving boundary modeling to account explicitly for the mechanics of
hyperelastic rods to capture a directional swelling of the underlying
elastomer.

Keywords: Moving boundary problem, finite element method, method of lines, a
priori error estimate, a posteriori error estimate, diffusion of chemicals
into rubber.

Mathematics Subject Classifications (2020). 65M15, 65M20, 65M60, 35R37

##  1 Introduction

Sharp interfaces moving in an a priori unknown way inside materials play a key
role in a number of study cases in science and technology, including in the
forecast of the durability of cementitious-based materials (cf. e.g. [8, 24,
25, 35]), large-time behavior of chemical species from the environment slowly
penetrating by diffusion and swelling rubber-based materials (cf. e.g. [2, 16,
26]), to controlling phase transitions like melting and freezing or solid-
solid changes in concrete (cf. e.g. [4, 29, 30]), to mention but a few. Due to
the inherent non-linearity of such moving boundary problems, analytical
representations of solutions are often either unavailable or not computable.
Hence, one has to rely on direct computational approaches to get insight for
instance in the behavior of large times of such moving sharp interfaces, as
this usually defines the lifetime of the material under investigation.

In the framework of this paper, we study a semi-discrete finite element
approximation of weak solutions to a one dimensional moving boundary problem
that models the diffusion of solvent into rubber (see Section 2). This is a
follow-up study of our recent work [26], where we proposed a finite element
approximation of solutions to a moving boundary problem which we used to
recover experimental data. Now, we explore the quality of our approximation
scheme. Specifically, we report on both a priori and a posteriori error
estimates for the mass concentration of the diffusants, and respectively, for
the position of the moving boundary. Our working techniques include integral
and energy-based estimates for the corresponding nonlinear parabolic problem
posed in a transformed fixed domain, combined with a suitable use of the
interpolation-trace inequality to handle the interface terms. At the technical
level, we were very much inspired by the references: [7, 11, 27, 28], and
[24]. It is worth noting that similar work has been done in related contexts.
For instance, in [8], the authors show the convergence of a numerical scheme
obtained by combining an Euler discretization in time with a Scharfetter-
Gummel discretization in space for a concrete carbonation model with moving
boundary reformulated for a fixed space domain. In [35], A. Zurek studies the
long time regime of the moving interface driving the concrete carbonation
reaction model by tailoring an implicit in time and finite volume in space
scheme. He proves that the approximate free boundary increases in time with
$\sqrt{t}$-law as theoretically predicted in [3]. In [22], one develops an
adaptive moving mesh method for the numerical solution of an enthalpy
formulation of a class of heat-conduction problems with phase change. The main
aim of [12] is to provide a comparison of several numerical methods including
displacing level sets, moving grids, and diffusing phase fields to address two
well-known Stefan problems arising as best formulations for phase
transformations like melting of a pure phase and diffusional solid-state phase
changes in binary systems.

To handle our problem, we decided to use the finite element method as this
fits best to the regularity of the (weak) solutions to our moving boundary
problem. Mind though that other discretization methods are likely to be
applicable as well. As our work is purely in 1D and no expensive computations
are expected, and as, on top of this, we wish to rely on open source
facilities, we chose Python for the implementation work.

We present here a preliminary investigation of this class of problems. This is
necessary before extending the current moving boundary modeling to account
explicitly for the mechanics of hyperelastic rods to capture a directional
swelling of the underlying elastomer. In this spirit, a natural next step
would be to perform the numerical analysis of a two-scale finite element
approximation of the setup described in [2].

The outline of this study is as follows: We formulate our moving boundary
problem in Section 2. The discussion of the setting of the model equations is
based on [26]. We collect in Section 3 our basic assumptions on parameters and
model components, as well as notations and existing preliminary results.
Section 4 contains the fixed domain transformation of our problem and the
definition of our concept of weak solutions which is then the subject of error
approximation estimates investigated here. Benefiting of the mathematical
analysis done for our problem in [17, 18], we are able to prove the global
existence of weak solutions to the semi-discrete problem and obtain the needed
uniform boundedness results to produce convergent numerical schemes. As main
result, we obtain a priori and a posteriori error estimates as listed in
Section 5. A couple of numerical experiments are discussed in Section 6.
Essentially, they support numerically the available experimental results.
Finally, a brief conclusion of this work is outlined in Section 7.

##  2 Model equations

We consider a thin slab of a dense rubber, denoted by $\Omega$ of vertical
length $L>0$, placed in contact with a diffusant reservoir. When the diffusant
concentration at the bottom face of the rubber exceeds some threshold, the
diffusant moves into the rubber creating a sharp interface that separates the
rubber $\Omega$ into two parts, the diffusant free region and diffusant-
penetrated region. Our region of interest is the diffusant-penetrated part
where the diffusant’s flux is assumed to satisfy Fick’s law. The actual
problem is to find the diffusant concentration profile inside the diffusant-
penetrated region and the location of the moving interface separating the
penetrated from the not-yet penetrated region. Such a setting is referred to
as a one-phase moving boundary problem. Formulations as a two-phase boundary
problem are possible as well, but are currently not in our focus; see e.g.
[10] for a nicely written  textbook regarding modeling with moving interfaces.

In this work, the modeling domain is the one–dimensional slab shown in Figure
1, which is the longitudinal line where $0<s(0)\leq s(t)\leq L$.

$0$$s(0)$$s(T_{f})$$L$Zone inside rubber with penetration of diffusants at
$t=0$Diffusant-free rubber at $t=0$ Figure 1: Sketch of one dimensional
geometry – a macroscopic thin slab made of rubber.

For a fixed observation time $T_{f}\in(0,\infty)$, the interval $[0,T_{f}]$ is
the time span of the process we are considering. Let $x\in[0,s(t)]$ and
$t\in[0,T_{f}]$ denote the space and respectively time variable, and let
$m(t,x)$ be the concentration of diffusant placed in position $x$ at time $t$.
The diffusants concentration $m(t,x)$ acts in the region $Q_{s}(T_{f})$
defined by

$Q_{s}(T_{f}):=\\{(t,x)|t\in(0,T_{f})\;\text{and}\;x\in(0,s(t))\\}.$

The problem reads: Find $m(t,x)$ and the position of the moving interface
$x=s(t)$ for $t\in(0,T_{f})$ such that the couple $(m(t,x),s(t))$ satisfies
the following

$\displaystyle\displaystyle\frac{\partial m}{\partial
t}-D\frac{\partial^{2}m}{\partial x^{2}}=0\;\;\;\
\text{in}\;\;\;Q_{s}(T_{f}),$ (1) $\displaystyle-D\frac{\partial m}{\partial
x}(t,0)=\beta(b(t)-\text{H}m(t,0))\;\;\;\text{for}\;\;t\in(0,T_{f}),$ (2)
$\displaystyle-D\frac{\partial m}{\partial
x}(t,s(t))=s^{\prime}(t)m(t,s(t))\;\;\;\text{for}\;\;t\in(0,T_{f}),$ (3)
$\displaystyle s^{\prime}(t)=a_{0}(m(t,s(t))-\sigma(s(t))\;\;\;\;\text{for
}\;\;\;t\in(0,T_{f}),$ (4) $\displaystyle
m(0,x)=m_{0}(x)\;\;\;\text{for}\;\;\;x\in[0,s(0)],$ (5) $\displaystyle
s(0)=s_{0}>0\;\text{with}\;\;0<s_{0}<s(t)<L,$ (6)

where $a_{0}>0$ is a kinetic coefficient, $\beta$ is a positive constant,
$D>0$ is a diffusion constant, $\text{H}>0$ is the Henry’s constant, $\sigma$
is a function on $\mathbb{R}$, $b$ is a given boundary function on $[0,T]$,
and $s_{0}>0$ is the initial position of the free boundary and $m_{0}$ is the
initial concentration of the diffusant.

The boundary condition (3) describes the mass conservation of diffusant
concentration at the moving boundary. It indicates that the diffusion
mechanism is responsible for pushing the interface. In particular (4) points
out that the mechanical behaviour (here it is about the swelling of the
rubber) also contributes to the motion of the moving penetration front. The
explanation of the model equations and the physical meaning of the parameters
are given in [26].

##  3 Notations, assumptions and preliminaries

In this section, we list our basic assumptions on the data, notations as well
as approximation properties of functions that are required for the error
analysis discussed in the next sections.

###  3.1 Function spaces and elementary inequalities

Let $u,v:\Omega\to\mathbb{R}$ denote two generic functions. Let
$W^{r,p}(\Omega)$ be the Sobolev space on domain $\Omega$ for $1\leq
p\leq\infty$ and $r\geq 0$. For $r=0$, we simply write $L^{p}(\Omega)$ in
place of $W^{0,p}(\Omega)$ with the norm $\|\cdot\|_{L^{p}(\Omega)}$ defined
as follows:

$\|u\|_{L^{p}(\Omega)}:=\begin{cases}\left(\displaystyle\int_{\Omega}|u(x)|^{p}dx\right)^{\frac{1}{p}}\;\;&\text{for}\;\;1\leq
p<\infty,\\\ \text{ess
sup}\\{|u(x)|:x\in\Omega\\}\;\;&\text{for}\;\;p=\infty,\end{cases}$

For $p=2$ and $r\geq 1$, we write $H^{r}(\Omega)$ in place of
$W^{r,2}(\Omega)$ with the norm $\|\cdot\|_{H^{r}(\Omega)}$ defined by

$\displaystyle\|u\|_{H^{r}(\Omega)}=\left(\sum_{|\alpha|\leq
r}\int_{\Omega}|\partial^{\alpha}u|^{2}dx\right)^{\frac{1}{2}}.$ (7)

In (7) $\partial^{\alpha}u$ denotes the $\alpha$’th derivative of $u$ in the
weak sense. Furthermore, for $L^{2}(\Omega)$ and $H^{r}(\Omega)$ we have the
following inner products.

$\displaystyle(u,v)_{L^{2}(\Omega)}:=\int_{\Omega}u(x)v(x)dx,$
$\displaystyle(u,v)_{H^{r}(\Omega)}:=\sum_{|\alpha|\leq
r}(\partial^{r}u,\partial^{r}v)_{L^{2}(\Omega)}.$

Let $X$ be a Banach space with norm $\|\cdot\|_{X}$ and $v:[0,T]\rightarrow X$
be a function. Correspondingly, $L^{p}(0,T,X)$ is a Bochner space endowed with
the norms

$\displaystyle\|v\|_{L^{p}(0,T,X)}:=\begin{cases}\left(\displaystyle\int_{0}^{T}\|v(\tau)\|_{X}^{p}d\tau\right)^{\frac{1}{p}}\;\;\;\;\;&\text{for}\;\;\;1\leq
p<\infty,\\\ \displaystyle\sup_{0\leq\tau\leq
T}\|v(\tau)\|_{X}\;\;\;&\text{for}\;\;p=\infty.\end{cases}$

More information on Sobolev and Bochner spaces with their various norms and
inner products can be found for instance in [1, 15]. For the convenience of
writing, we denote $u(t,0)$ and $u(t,1)$ by $u(0)$ and $u(1)$, respectively.
We also use the prime $(^{\prime})$ to point out the derivative with respect
to time variable, and $\|\cdot\|$ and $(\cdot,\cdot)$ for the norm and,
respectively, inner product in $L^{2}(\Omega)$. Furthermore,
$\|\cdot\|_{\infty}$ refers to the norm of $L^{\infty}(\Omega)$.
We list a few elementary inequalities that we frequently use in this work.

  1. (i)

Young’s inequality:

$\displaystyle ab\leq\xi a^{p}+c_{\xi}b^{q},$ (8)

where
$a,b\in\mathbb{R}_{+},\;\xi>0,\;c_{\xi}:=\displaystyle\frac{1}{q}\frac{1}{\sqrt[p]{(\xi
p)^{q}}}>0,\;\displaystyle\frac{1}{p}+\frac{1}{q}=1$ and $p\in(1,\infty).$

  2. (ii)

Interpolation inequality: For all $u\in H^{1}(0,1)$, there exists a constant
$\hat{c}>0$ depending on $\theta\in[\frac{1}{2},1)$ such that

$\displaystyle\|u\|_{\infty}\leq\hat{c}\|u\|^{\theta}\left\|u\right\|_{H^{1}(0,1)}^{1-\theta}.$
(9)

For $\theta=1/2$, one gets

$\displaystyle\|u\|_{\infty}^{2}\leq\hat{c}\left(\xi\left\|\frac{\partial
u}{\partial y}\right\|^{2}+(\xi+c_{\xi})\|u\|^{2}\right),$

where $\xi$ and $c_{\xi}$ are as in (8). See details in [34] p. 285 (example
21.62).

###  3.2 Assumptions on parameters

Throughout this paper, we assume the following restrictions on the parameters.

  1. (A1)

$a_{0},\;\text{H},\;D,\;s_{0},\;T_{f}$ are positive constants.

  2. (A2)

$b\in W^{1,2}(0,T_{f})$ with $0<b_{*}\leq b\leq b^{*}$ on $(0,T_{f})$, where
$b_{*}$ and $b^{*}$ are positive constants.

  3. (A3)

$\beta\in C^{1}(\mathbb{R})\cap W^{1,\infty}(\mathbb{R})$ such that $\beta=0$
on $(\infty,0]$, and there exists $r_{\beta}>0$ such that $\beta^{\prime}>0$
on $(0,r_{\beta})$ and $\beta=k_{0}$ on $[r_{\beta},+\infty)$, where
$k_{0}>0$.

  4. (A4)

$\sigma\in C^{1}(\mathbb{R})\cap W^{1,\infty}(\mathbb{R})$ such that
$\sigma=0$ on
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(-\infty,0)}$,
and there exists $r_{\sigma}$ such that $\sigma^{\prime}>0$ on
$(0,r_{\sigma})$ and $\sigma=c_{0}$ on $[r_{\sigma},+\infty),$ where $c_{0}$
satisfies

$\displaystyle 0<c_{0}<\min\\{2\sigma(0),b^{*}\text{H}^{-1}\\}.$ (10)

  5. (A5)

$0<s_{0}<r_{\sigma}$ and $m_{0}\in H^{1}(0,s_{0})$ such that $\sigma(0)\leq
u_{0}\leq b^{*}\text{H}^{-1}$ on $[0,s_{0}].$

The assumptions (A1)–(A5) are adopted from [18], where the authors have proved
the global solvability of the problem and continuous dependence estimates of
the solution with respect to the initial data.

###  3.3 Basic facts from approximation theory

Let $N\in\mathbb{N}$ be given. We set $0=y_{0}<y_{1}<\cdots<y_{N-1}=1$ as
discretization points in the interval $[0,1]$. We set $k_{i}:=y_{i+1}-y_{i}$
for $i\in\\{0,1,\cdots,N-2\\}$ and $k:=\displaystyle\max\\{k_{i}:\
i\in\\{0,1,\cdots,N-2\\}\\}$. We introduce the space

$\displaystyle V_{k}:=\\{\nu\in
C[0,1]:\nu|_{[y_{j},y_{j+1}]}\in\mathbb{P}_{1}\\},$ (11)

where $\mathbb{P}_{1}$ represents the set of polynomials of degree one. Let
$\\{\phi_{i}\\}_{i=0}^{N-1}$ be the set of basis functions for the space
$V_{k}$ defined by

$\phi_{i}(y)=\begin{cases}0&\;\;\;\;\;\;\text{if}\;\;\;\;y<y_{i-1}\\\
\displaystyle\frac{y-y_{i-1}}{k_{i-1}}&\;\;\;\;\;\;\text{if}\;\;\;y_{i-1}\leq
y<y_{i}\\\
\displaystyle\frac{y_{i+1}-y}{k_{i}}&\;\;\;\;\;\;\text{if}\;\;\;y_{i}\leq
y<y_{i+1}\\\ 0&\;\;\;\;\;\;\text{if}\;\;\;\;y_{i+1}\leq y.\\\ \end{cases}$

We define the interpolation operator $I_{k}:C[0,1]\rightarrow V_{k}$ by

$\displaystyle(I_{k}u)(y):=\sum_{i=0}^{N-1}u(y_{i},t)\phi_{i}(y).$

Here the function $I_{k}u$ is called the Lagrange interpolant of $u$ of degree
1; for more details see e.g. [20].

######  Lemma 3.1.

Take $\theta\in[\frac{1}{2},1)$ and $\psi\in H^{2}(0,1)$. Then there exist
strictly positive constants $\gamma_{1},\;\gamma_{2}$ and $\gamma_{3}$ such
that the Lagrange interpolant $I_{k}\psi$ of $\psi$ satisfies the following
estimates:

  1. (i)

$\|\psi-I_{k}\psi\|\leq\gamma_{1}k^{2}\|\psi\|_{H^{2}(0,1)}$

  2. (ii)

$\left\|\displaystyle\frac{\partial}{\partial y}(\psi-
I_{k}\psi)\right\|\leq\gamma_{2}k\|\psi\|_{H^{2}(0,1)}$

  3. (iii)

$|\psi(0)-I_{k}\psi(0)|\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{c}\left(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta}\right)\|\psi\|_{H^{2}(0,1)}}$

  4. (iv)

$|\psi(1)-I_{k}\psi(1)|\leq\hat{c}\left(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta}\right)\|\psi\|_{H^{2}(0,1)}$

###### Proof.

The inequalities (i) and (ii) are standard results. For details on their
proof, see for instance page 61 in [20] and page 3 in [32]. To show (iii), we
use the interpolation inequality (9) together with (i) and (ii) to obtain

$\displaystyle|\psi(0)-I_{k}\psi(0)|$ $\displaystyle\leq\hat{c}\|\psi-
I_{k}\psi\|^{\theta}_{L^{2}(0,1)}\left\|\psi-
I_{k}\psi\right\|^{1-\theta}_{H^{1}(0,1)}$ $\displaystyle\leq\hat{c}\|\psi-
I_{k}\psi\|^{\theta}_{L^{2}(0,1)}\left(\|\psi-
I_{k}\psi\|^{1-\theta}+\left\|\displaystyle\frac{\partial}{\partial y}(\psi-
I_{k}\psi)\right\|^{1-\theta}\right)$
$\displaystyle\leq\hat{c}\left(\gamma_{1}k^{2}+\gamma_{1}^{\theta}\gamma_{2}^{1-\theta}k^{1+\theta}\right)\|\psi\|_{H^{2}(0,1)}.$

Taking $\gamma_{3}:=\gamma_{1}^{\theta}\gamma_{2}^{1-\theta}$ gives the
estimate $({\rm\ref{l3}})$. A similar argument applied to $\psi(1)$ gives
$({\rm\ref{l4}})$. ∎

##  4 Fixed-domain transformation and definition of weak solutions

Firstly, we perform the non-dimensionalization of the model equations (1)–(6).
We then transform the non-dimensional model equations from the a priori
unknown non-cylinderical domain into the cylinderical domain
$Q(T):=\\{(\tau,y)|\;\tau\in(0,T)\;\text{and}\;y\in(0,1)\\}$ by using the
Landau transformation $y=x/s(t)$, see for instance [19]. For more details on
non-dimensionalization and transformation, we refer the reader to [26] where
the preliminary steps are done. In dimensionless form, the transformed problem
reads as follows:

$\displaystyle\displaystyle\frac{\partial{u}}{\partial{\tau}}-y\frac{h^{\prime}(\tau)}{h(\tau)}\frac{\partial{u}}{\partial{y}}-\frac{1}{(h(\tau))^{2}}\frac{\partial^{2}{u}}{\partial{y^{2}}}=0\;\;\;\
\text{in}\;\;\;Q(T),$ (12) $\displaystyle-\frac{1}{h(\tau)}\frac{\partial
u}{\partial
y}(\tau,0)=\text{\rm{Bi}}\left(\frac{b(\tau)}{m_{0}}-\text{H}u(\tau,0)\right)\;\;\;\text{for}\;\;\tau\in(0,T),$
(13) $\displaystyle-\frac{1}{h(\tau)}\frac{\partial u}{\partial
y}(\tau,1)=h^{\prime}(\tau)u(\tau,1)\;\;\;\text{for}\;\;\tau\in(0,T),$ (14)
$\displaystyle
h^{\prime}(\tau)=A_{0}\left(u(\tau,1)-\frac{\sigma(h(\tau))}{m_{0}}\right)\;\;\;\text{for}\;\;\tau\in(0,T)$
(15) $\displaystyle u(0,y)=u_{0}(y)\;\;\;\text{for}\;\;y\in[0,1],$ (16)
$\displaystyle h(0)=h_{0}.$ (17)

We refer to the system (12)–(17) posed in the cylinderical domain $Q(T)$ as
problem $(P)$.

######  Remark 4.1.

We refer the reader to [26] for the definition of dimensionless quantities
$u,\;h,\;\tau,\;y,\;T,\;{\rm Bi},\;A_{0}$. Here we only mention that Bi is the
mass transfer Biot number and $A_{0}$ is the Thiele modulus.

######  Definition 4.1.

(Weak Solution to ($P$)). We call the couple $(u,h)$ a weak solution to
problem ($P$) on $S_{T}:=(0,T)$ if and only if

$\displaystyle h\in W^{1,\infty}(S_{T})\;\;\text{with}\;\;h_{0}<h(T)\leq L,$
$\displaystyle u\in W^{1,2}(Q(T))\cap L^{\infty}(S_{T},H^{1}(0,1))\cap
L^{2}(S_{T},H^{2}(0,1)),$

such that for all $\tau\in S_{T}$ the following relations hold

$\displaystyle\displaystyle\left(\frac{\partial
u}{\partial\tau},\varphi\right)$
$\displaystyle-\frac{h^{\prime}(\tau)}{h(\tau)}\left(y\frac{\partial
u}{\partial
y},\varphi\right)+\frac{1}{(h(\tau))^{2}}\left(\frac{\partial{u}}{\partial{y}},\frac{\partial\varphi}{\partial
y}\right)$ $\displaystyle-\frac{1}{h(\tau)}{\rm
Bi}\left(\frac{b(\tau)}{m_{0}}-{\rm
H}u(\tau,0)\right)\varphi(0)+\frac{h^{\prime}(\tau)}{h(\tau)}u(\tau,1)\varphi(1)=0\;\;\text{for
all}\;\varphi\in H^{1}(0,1),$ (18) $\displaystyle
h^{\prime}(\tau)=A_{0}\left(u(\tau,1)-\frac{\sigma(h(\tau))}{m_{0}}\right),$
(19) $\displaystyle u(0,y)=u_{0}(y)\;\;for\;\;y\in[0,1],$ (20) $\displaystyle
h(0)=h_{0}.$ (21)

######  Theorem 4.1.

If (A1)–(A5) hold, then problem $(P)$ has a unique solution $(u,h)$ on $S_{T}$
in the sense of Definition 4.1.

###### Proof.

We refer the reader to Theorem 2.4 in [17] for a statement of the local
existence of weak solutions to problem $(P)$ and to Theorem 3.3 and Theorem
3.4 in [18] for a way to ensure the global existence and continuous dependence
with respect to initial data.∎

We now define the finite element Galerkin approximation to (18)–(21) on the
finite dimensional subspace $V_{k}$. The semi-discrete approximation $u_{k}$
and $h_{k}$ of $u$ and $h$ is now defined to be the mapping
$u_{k}:[0,T]\rightarrow V_{k}$ and $h_{k}:[0,T]\rightarrow\mathbb{R}_{+}$ such
that (22)–(25) holds. We denote the semi-discrete form (22)–(25) of problem
$(P)$ by $(P_{d})$.

######  Definition 4.2.

(Weak Solution to $(P_{d})$). We call the couple $(u_{k},h_{k})$ a weak
solution to problem ($P_{d}$) if and only if there is a $S_{{T}}:=(0,{T})$
(for some ${T}>0$) such that

$\displaystyle h_{k}\in
W^{1,\infty}(S_{{T}})\;\;\text{with}\;\;h_{0}<h_{k}({T})\leq L$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}u_{k}\in
H^{1}(S_{{T}},V_{k})\cap L^{2}(S_{{T}},H^{1}(0,1))\cap
L^{\infty}(S_{{T}},L^{2}(0,1))}$

and for all $\tau\in S_{{T}}$ it holds

$\displaystyle\displaystyle\left(\frac{\partial
u_{k}}{\partial\tau},\varphi_{k}\right)$
$\displaystyle-\frac{h_{k}^{\prime}(\tau)}{h_{k}(\tau)}\left(y\frac{\partial
u_{k}}{\partial
y},\varphi_{k}\right)+\frac{1}{(h_{k}(\tau))^{2}}\left(\frac{\partial{u_{k}}}{\partial{y}},\frac{\partial\varphi_{k}}{\partial
y}\right)$ $\displaystyle-\frac{1}{h_{k}(\tau)}{\rm
Bi}\left(\frac{b(\tau)}{m_{0}}-{\rm
H}u_{k}(\tau,0)\right)\varphi_{k}(0)+\frac{h_{k}^{\prime}(\tau)}{h_{k}(\tau)}u_{k}(\tau,1)\varphi_{k}(1)=0\;\;\text{for
all}\;\;\varphi_{k}\in V_{k},$ (22) $\displaystyle
h_{k}^{\prime}(\tau)=A_{0}\left(u_{k}(\tau,1)-\frac{\sigma(h_{k}(\tau))}{m_{0}}\right),$
(23) $\displaystyle u_{k}(0)=u_{0,k}(y)\;\;{\rm for}\;\;y\in[0,1],$ (24)
$\displaystyle h_{k}(0)=h_{0}.$ (25)

######  Lemma 4.1.

Assume (A1)–(A5) hold. Then there exist a time $\hat{T}\in(0,T]$ and positive
constants $L,M_{1},M_{2}$ (not depending of $k$) such that for a.e.
$\tau\in(0,\bar{T})$ the following inequalities hold true for the pair
$(u_{k},h_{k})$ arising in Definition 4.2:

  1. (i)

$0<h_{0}\leq h_{k}(\tau)\leq L$

  2. (ii)

$0<u_{k}(\tau,y)<M_{1}$

  3. (iii)

$|h_{k}^{\prime}(\tau)|\leq M_{2}$.

###### Proof.

(i) is built in the concept of weak solution detailed in Definition 4.2. It
does not require a proof. We added it here simply to stress the importance of
the fact the we work exclusively in a bounded moving domain. (ii) is the main
statement here. This holds true as a consequence of the fact that the space
continuous version of the statement (i.e. $0<u(\tau,y)<M_{1}$) holds true; we
rely on the arguments of the proof of Theorem 3.1 in [18], combined with the
fact that the treated geometry is one dimensional. Hence, $\hat{T}>0$ is
possibly small, which is sufficient for deriving our next results. Note though
that a discrete version of the Stampacchia trick, worked out with details in
[14], can potentially be applied here as well in order to replace the local
time $\hat{T}$ with a maximal time. Alternative arguments employing the
structure of the problem as in [33] or based on linear simplicial finite
elements as in [6] can also be used in principle. (iii) is a direct
consequence of (i) and (ii) combined with (23).

∎

######  Theorem 4.2.

Let the hypothesis of Lemma 4.1 be fulfilled. Then it exists a unique solution

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(u_{k},h_{k})\in
H^{1}(S_{\hat{T}},V_{k})\cap L^{2}(S_{\hat{T}},H^{1}(0,1))\cap
L^{\infty}(S_{\hat{T}},L^{2}(0,1))\times W^{1,\infty}(S_{\hat{T}})}$

in the sense of Definition 4.2. Furthermore, there exists a constant
$\tilde{c}>0$ (independent of $k$) such that

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\max_{0\leq\tau\leq\hat{T}}}\|u_{k}\|_{L^{2}(0,1)}^{2}+\int_{0}^{\hat{T}}\left\|\frac{\partial
u_{k}}{\partial y}\right\|^{2}_{L^{2}(0,1)}d\tau\leq\tilde{c}.$ (26)

###### Proof.

Let $V_{k}$ be the finite dimensional subspace defined in (11) constructed
based on the span of the hat functions
$\\{\phi_{j}\\},j\in\\{0,1,\cdots,N-1\\}$. Let
$\alpha_{j}:(0,\hat{T})\rightarrow\mathbb{R}$ denote the Galerkin projection
coefficient for $j$th degree of freedom. Then the finite-dimensional Galerkin
approximation of the function $u$ is defined by

$\displaystyle u_{k}(\tau,y):=\sum_{j=0}^{N-1}\alpha_{j}(\tau)\phi_{j}(y),$

where the coefficients $\alpha_{j}(\tau),\;j\in\\{0,1,\dots N-1\\}$ are
determined by the following relations:

$\displaystyle\displaystyle\left(\frac{\partial
u_{k}}{\partial\tau},\varphi_{k}\right)-\frac{h_{k}^{\prime}(\tau)}{h_{k}(\tau)}\left(y\frac{\partial
u_{k}}{\partial y},\varphi_{k}\right)$
$\displaystyle+\frac{1}{(h_{k}(\tau))^{2}}\left(\frac{\partial{u_{k}}}{\partial{y}},\frac{\partial\varphi_{k}}{\partial
y}\right)$ $\displaystyle-\frac{1}{h_{k}(\tau)}{\rm
Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u_{k}(\tau,0)\right)\varphi_{k}(0)+\frac{h_{k}^{\prime}(\tau)}{h_{k}(\tau)}u_{k}(\tau,1)\varphi_{k}(1)=0,$
(27) $\displaystyle
h_{k}^{\prime}(\tau)=A_{0}\left(u_{k}(\tau,1)-\frac{\sigma(h_{k}(\tau))}{m_{0}}\right),\;\;\;\tau\in(0,\hat{T})$
(28)

for all $\varphi_{k}\in\text{span}\\{\phi_{j}\\},\;j\in\\{0,1,\cdots,N-1\\}$
and

$\displaystyle\alpha_{j}(0)=(u_{0,k},\phi_{j}),$ (29) $\displaystyle
h_{k}(0)=h_{0}.$ (30)

Taking in (27) and (28) as test function $\varphi_{k}=\phi_{j}$ for
$j\in\\{0,1,\cdots,N-1\\}$, we obtain the following system of ordinary
differential equations for the unknown
$\alpha=(\alpha_{j})_{j=0,1,\cdots,N-1}$ and $h_{k}$:

$\displaystyle\sum_{i=0}^{N-1}M_{i}\alpha_{i}^{\prime}(\tau)-\frac{h_{k}^{\prime}}{h_{k}}\sum_{i=0}^{N-1}K_{i}\alpha_{i}+\frac{1}{h_{k}^{2}}\sum_{i=0}^{N-1}A_{i}\alpha_{i}=\frac{1}{h_{k}}{\rm
Bi}\left(\frac{b(\tau)}{m_{0}}\phi(0)-\text{H}\alpha\right)-\frac{h_{k}^{\prime}}{h_{k}}\alpha=:G_{1}(\alpha,h_{k}),$
(31) $\displaystyle
h_{k}^{\prime}(\tau)=A_{0}\left(\sum_{i=0}^{N-1}\alpha_{i}\phi_{i}(1)-\frac{\sigma(h_{k}(\tau))}{m_{0}}\right)=:G_{2}(\alpha,h_{k}),$
(32)

where

$\displaystyle(M_{i})_{j}:=\int_{0}^{1}\phi_{i}\phi_{j}dy,$ (33)
$\displaystyle(K_{i})_{j}:=\int_{0}^{1}y\frac{\partial\phi_{i}}{\partial
y}\phi_{j}dy,$ (34)
$\displaystyle(A_{i})_{j}:=\int_{0}^{1}\frac{\partial\phi_{i}}{\partial
y}\frac{\partial\phi_{j}}{\partial y}dy.$ (35)

Firstly, we prove that $G_{2}$ is Lipschitz. Let $(\alpha,h_{k})$ and
$(\beta,\tilde{h}_{k})$ be two pairs.

$\displaystyle\left|G_{2}(\alpha,h_{k})-G_{2}(\beta,\tilde{h}_{k})\right|\leq
A_{0}\left(\sum_{i=0}^{N-1}\left|\alpha_{i}(\tau)-\beta_{i}(\tau)\right|\left|\phi_{i}(1)\right|+\frac{1}{m_{0}}\left|\sigma(h_{k}(\tau))-\sigma(\tilde{h}_{k}(\tau))\right|\right).$
(36)

Using (A4) in (36), we get

$\displaystyle\left|G_{2}(\alpha,h_{k})-G_{2}(\beta,\tilde{h}_{k})\right|$
$\displaystyle\leq
A_{0}\left(\sum_{i=0}^{N-1}\left|\alpha_{i}(\tau)-\beta_{i}(\tau)\right|\left|\phi_{i}(1)\right|+\frac{\mathcal{L}}{m_{0}}\left|h_{k}(\tau)-\tilde{h}_{k}(\tau)\right|\right)$
$\displaystyle\leq\mathcal{M}\left(\sum_{i=0}^{N-1}\left|\alpha_{i}(\tau)-\beta_{i}(\tau)\right|+\left|h_{k}(\tau)-\tilde{h}_{k}(\tau)\right|\right)$
$\displaystyle=\mathcal{M}|(\alpha,h_{k})-(\beta,\tilde{h}_{k})|,$

where $\mathcal{L}$ is a Lipschitz constant and

$\mathcal{M}:=\max\left\\{A_{0}\max_{0\leq i\leq
N-1}|\phi_{i}(1)|,\frac{A_{0}\mathcal{L}}{m_{0}}\right\\}.$

Thus, $G_{2}$ is Lipschitz. Now, we show that $G_{1}$ is Lipschitz.

$\displaystyle G_{1}(\alpha,h_{k})-G_{1}(\beta,\tilde{h}_{k})={\rm
Bi}\frac{b(\tau)}{m_{0}}\left(\frac{1}{h_{k}}-\frac{1}{\tilde{h}_{k}}\right)\phi(0)-{\rm
Bi}\,\text{H}\left(\frac{\alpha}{h_{k}}-\frac{\beta}{\tilde{h}_{k}}\right)-\left(\frac{h_{k}^{\prime}}{h_{k}}\alpha-\frac{\tilde{h}_{k}^{\prime}}{\tilde{h}_{k}}\beta\right).$
(37)

Using (A2) in (37) yields

$\displaystyle\left|G_{1}(\alpha,h_{k})-G_{1}(\beta,\tilde{h}_{k})\right|\leq$
$\displaystyle\;\text{Bi}\frac{b^{*}}{m_{0}h_{k}\tilde{h}_{k}}|h_{k}-\tilde{h}_{k}||\phi(0)|+\text{Bi}\,\text{H}\left|\frac{\alpha}{h_{k}}-\frac{\beta}{\tilde{h}_{k}}\right|$
$\displaystyle+\left|\frac{h_{k}^{\prime}}{h_{k}}\alpha-\frac{\tilde{h}_{k}^{\prime}}{\tilde{h}_{k}}\beta\right|$
$\displaystyle=$ $\displaystyle\sum_{\ell=1}^{3}I_{\ell},$

where

$\displaystyle I_{1}$
$\displaystyle:=\text{Bi}\frac{b^{*}}{m_{0}h_{k}\tilde{h}_{k}}|h_{k}-\tilde{h}_{k}||\phi(0)|\leq\text{Bi}\frac{b^{*}}{m_{0}h_{k}\tilde{h}_{k}}|h_{k}-\tilde{h}_{k}|,$
$\displaystyle I_{2}$
$\displaystyle:=\text{Bi}\,\text{H}\left|\frac{\alpha}{h_{k}}-\frac{\beta}{\tilde{h}_{k}}\right|\leq\text{Bi}\,\text{H}\left(|\alpha|\frac{|h_{k}-\tilde{h}_{k}|}{h_{k}\tilde{h}_{k}}+\frac{|\alpha-\beta|}{\tilde{h}_{k}}\right),$
$\displaystyle I_{3}$
$\displaystyle:=\left|\frac{h_{k}^{\prime}}{h_{k}}\alpha-\frac{\tilde{h}_{k}^{\prime}}{\tilde{h}_{k}}\beta\right|$
$\displaystyle=\left|h_{k}^{\prime}\left(\frac{\alpha}{h_{k}}-\frac{\beta}{\tilde{h}_{k}}\right)+\frac{\beta}{\tilde{h}_{k}}\left(h_{k}^{\prime}-\tilde{h}_{k}^{\prime}\right)\right|$
$\displaystyle\leq\left|h_{k}^{\prime}\right|\left|\frac{\alpha}{h_{k}}-\frac{\beta}{\tilde{h}_{k}}\right|+\mathcal{L}\frac{|\beta|}{|\tilde{h}_{k}|}\left|h_{k}-\tilde{h}_{k}\right|$
$\displaystyle\leq\left|h_{k}^{\prime}\right|\left(|\alpha|\frac{|h_{k}-\tilde{h}_{k}|}{h_{k}\tilde{h}_{k}}+\frac{|\alpha-\beta|}{\tilde{h}_{k}}\right)+\mathcal{L}\frac{|\beta|}{|\tilde{h}_{k}|}\left|h_{k}-\tilde{h}_{k}\right|.$

This shows that $G_{1}$ is Lipschitz continuous. By standard arguments for
systems of ordinary differential equations, the problem (29)-(32) has a unique
solution

$(\alpha,h_{k})\in C^{1}([0,\hat{T}])^{N}\times W^{1,\infty}(0,\hat{T}).$

We now prove the uniform estimate for the solution $u$ to the finite
dimensional problem.
Taking $\varphi=u_{k}$ in (27) yields

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\|u_{k}(\tau)\|^{2}$
$\displaystyle+\frac{1}{(h_{k}(\tau))^{2}}\left\|\frac{\partial
u_{k}(\tau)}{\partial
y}\right\|^{2}=\int_{0}^{1}\frac{h_{k}^{\prime}}{h_{k}}y\frac{\partial
u_{k}}{\partial y}u_{k}dy$
$\displaystyle+\frac{1}{h_{k}(\tau)}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u_{k}(\tau,0)\right)u_{k}(\tau,0)+\frac{h_{k}^{\prime}(\tau)}{h_{k}(\tau)}u_{k}(\tau,1)u_{k}(\tau,1).$
(38)

Using Hölder’s inequality for the first term on the right hand side of (38),
it holds that

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\|u_{k}(\tau)\|^{2}+\frac{1}{(h_{k}(\tau))^{2}}\left\|\frac{\partial
u_{k}(\tau)}{\partial y}\right\|^{2}$
$\displaystyle\leq\frac{|h_{k}^{\prime}|}{h_{k}}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|_{L^{2}(\Omega)}\left\|u_{k}\right\|_{L^{2}(\Omega)}$
$\displaystyle+\frac{\rm{Bi}}{h_{k}}\frac{b^{*}}{m_{0}}|u_{k}(\tau,0)|+\frac{\left|h_{k}^{\prime}\right|}{h_{k}}\left|u_{k}(\tau,1)\right|^{2}.$
(39)

We note here that, by the Sobolev’s embedding inequality in one space
dimension, it holds

$\displaystyle|\vartheta(\tau,y)|^{2}\leq
C_{e}\left\|\vartheta(\tau)\right\|_{H^{1}(0,1)}\left\|\vartheta(\tau)\right\|_{L^{2}(0,1)}\;\;\text{for}\;\vartheta\in
H^{1}(0,1)\;\text{and}\;y\in[0,1],$ (40)

where $C_{e}$ is a positive constant. Using (40), the third term on the right
hand side of (39) becomes

$\displaystyle\frac{\left|h_{k}^{\prime}\right|}{h_{k}}\left|u_{k}(\tau,1)\right|^{2}$
$\displaystyle\leq
C_{e}\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left\|u_{k}(\tau)\right\|_{H^{1}(0,1)}\left\|u_{k}(\tau)\right\|_{L^{2}(0,1)}$
$\displaystyle\leq
C_{e}\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left(\left\|\frac{\partial
u_{k}(\tau)}{\partial
y}\right\|_{L^{2}(0,1)}\left\|u_{k}(\tau)\right\|_{L^{2}(0,1)}+\left\|u_{k}(\tau)\right\|_{L^{2}(0,1)}^{2}\right).$
(41)

Using (41), (39) becomes

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\|u_{k}\|^{2}+\frac{1}{(h_{k})^{2}}\left\|\frac{\partial
u_{k}}{\partial y}\right\|^{2}\leq$
$\displaystyle(1+C_{e})\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left\|\frac{\partial
u_{k}}{\partial y}\right\|_{L^{2}(0,1)}\left\|u_{k}\right\|_{L^{2}(0,1)}$
$\displaystyle+C_{e}\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left\|u_{k}\right\|_{L^{2}(0,1)}^{2}+\frac{1}{h_{0}}\frac{b^{*}}{m_{0}}\|u_{k}\|_{H^{1}(0,1)}.$
(42)

Using Young’s inequality, (42) leads to

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\|u_{k}\|^{2}+\frac{1}{2L^{2}}\left\|\frac{\partial
u_{k}}{\partial y}\right\|^{2}\leq$
$\displaystyle(1+C_{e})\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left(\xi\left\|\frac{\partial
u_{k}}{\partial
y}\right\|_{L^{2}(0,1)}^{2}+c_{\xi}\left\|u_{k}\right\|_{L^{2}(0,1)}^{2}\right)$
$\displaystyle+C_{e}\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left\|u_{k}\right\|_{L^{2}(0,1)}^{2}+\xi\left\|\frac{\partial
u_{k}}{\partial y}\right\|_{L^{2}(0,1)}^{2}$
$\displaystyle+\xi\left\|u_{k}\right\|_{L^{2}(0,1)}^{2}+\frac{c_{\xi}}{h_{0}^{2}}\frac{(b^{*})^{2}}{m_{0}^{2}}.$

Finally, we get the following inequality

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\|u_{k}\|^{2}+M_{1}\left\|\frac{\partial
u_{k}}{\partial y}\right\|^{2}\leq
M_{2}\left\|u_{k}\right\|_{L^{2}(0,1)}^{2}+M_{3},$ (43)

where

$\displaystyle
M_{1}:=\frac{1}{2L^{2}}-\left((1+C_{e})\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}+1\right)\xi,$
$\displaystyle
M_{2}:=\frac{\left\|h_{k}^{\prime}\right\|_{L^{\infty}(S_{\hat{T}})}}{h_{0}}\left(c_{\xi}+C_{e}(c_{\xi}+1)\right)+\xi,$
$\displaystyle M_{3}:=\frac{c_{\xi}}{h_{0}^{2}}\frac{(b^{*})^{2}}{m_{0}^{2}}.$

Choosing a sufficiently small $\xi$ with $M_{1}>0$ and then applying
Gronwall’s inequality gives the following inequality holds

$\displaystyle\|u_{k}(\tau)\|^{2}\leq
c(\hat{T},h_{0},C_{e})\left(\|u_{k}(0)\|^{2}+M_{3}\hat{T}\right),$ (44)

for all $0\leq\tau\leq\hat{T}.$ Since $\|u_{k}(0)\|^{2}\leq\|u_{0,k}\|^{2}$,
(44) yields

$\displaystyle\max_{0\leq\tau\leq\hat{T}}\|u_{k}(\tau)\|^{2}\leq\tilde{c}.$
(45)

Integrating (43) from $0$ to $\hat{T}$ and employ the inequality (45) to get

$\displaystyle\int_{0}^{\hat{T}}\left\|\frac{\partial u_{k}}{\partial
y}\right\|^{2}d\tau\leq\tilde{c}.$

This concludes the proof of (26). ∎

######  Remark 4.2.

The entries of the matrices $M,\;K$ and $A$ given in (33), (34) and (35) are
computed explicitly benefiting of the structure of the basis elements
$\phi_{j}\in V_{k}$, usually piecewise polynomials of some preset degree
defined in $\Omega$; see [26] for the explicit form of the matrix $K$ and $A$
when using as basis piecewise linear functions.

##  5 Main results

In this Section, we prove a priori and a posteriori error estimates between
the weak solution to $(P)$ and weak solution to a semi-discrete version of
$(P)$. The discretization in space is done via the finite element method [20].

######  Theorem 5.1.

(A priori error estimate) Assume (A1)–(A5) hold. Additionally, take $u_{0}\in
H^{2}(0,1)$. Let $(u,h)$ and $(u_{k},h_{k})$ be the corresponding weak
solutions to problem $(P)$ and $(P_{d})$ in the sense of Definition 4.1 and
Definition 4.2, respectively. Then there exists a constant $c>0$ (not
depending on k) such that

$\displaystyle\|u-u_{k}\|^{2}_{L^{\infty}(S_{\hat{T}},L^{2}(0,1)))\cap
L^{2}(S_{\hat{T}},H^{1}(0,1))}+\|h-h_{k}\|^{2}_{H^{1}(S_{\hat{T}})}\leq
ck^{2}.$ (46)

###### Proof.

We consider the time interval $S_{\hat{T}}$ on which both continuous and
discrete solutions to (12)–(17) exist and are uniquely defined. Let
$e:=u-u_{k}$ and $h-h_{k}$ be the pointwise errors of the approximation. By
subtracting (22) from (18) and choosing $\varphi=v_{k}\in V_{k}$, we obtain
the following identity:

$\displaystyle\left(\frac{\partial u}{\partial\tau},v_{k}\right)-$
$\displaystyle\left(\frac{\partial
u_{k}}{\partial\tau},v_{k}\right)+\frac{1}{h^{2}}\left(\frac{\partial{u}}{\partial{y}},\frac{\partial
v_{k}}{\partial y}\right)-\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v_{k}}{\partial y}\right)$
$\displaystyle-\left(\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial
u}{\partial y}v_{k}dy-\frac{h_{k}^{\prime}}{h_{k}}\int_{0}^{1}y\frac{\partial
u_{k}}{\partial
y}v_{k}dy\right)+\frac{h^{\prime}}{h}u(\tau,1)v_{k}(1)-\frac{h_{k}^{\prime}}{h_{k}}u_{k}(\tau,1)v_{k}(1)$
$\displaystyle-\left(\frac{1}{h}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u(\tau,0)\right)v_{k}(0)-\frac{1}{h_{k}}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u_{k}(\tau,0)\right)v_{k}(0)\right)=0,$
(47)

which holds for all $v_{k}\in V_{k}$ and for almost every $\tau\in
S_{\hat{T}}$.
Arranging conveniently the terms in (47) yields

$\displaystyle\left(\frac{\partial e}{\partial\tau},v_{k}\right)$
$\displaystyle+\frac{1}{h^{2}}\left(\frac{\partial{e}}{\partial{y}},\frac{\partial
v_{k}}{\partial
y}\right)-\left(\frac{1}{h_{k}^{2}}-\frac{1}{h^{2}}\right)\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v_{k}}{\partial y}\right)$
$\displaystyle-\left(\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial
e}{\partial
y}v_{k}dy+\left(\frac{h^{\prime}}{h}-\frac{h_{k}^{\prime}}{h_{k}}\right)\int_{0}^{1}y\frac{\partial
u_{k}}{\partial y}v_{k}dy\right)$
$\displaystyle+\frac{h^{\prime}}{h}e(\tau,1)v_{k}(1)+\left(\frac{h^{\prime}}{h}-\frac{h_{k}^{\prime}}{h_{k}}\right)u_{k}(\tau,1)v_{k}(1)$
$\displaystyle-\left(\text{Bi}\frac{b(\tau)}{m_{0}}\left(\frac{1}{h}-\frac{1}{h_{k}}\right)v_{k}(0)-\text{Bi}\;\text{H}\left(\frac{u(\tau,0)}{h}-\frac{u_{k}(\tau,0)}{h_{k}}\right)v_{k}(0)\right)=0.$
(48)

In (48), we take as test function $v_{k}:=w_{k}-u_{k}\in V_{k}$ and use the
decomposition $v_{k}=(w_{k}-u)+e$. Then (48) becomes

$\displaystyle\left(\frac{\partial e}{\partial\tau},e\right)$
$\displaystyle+\left(\frac{\partial
e}{\partial\tau},w_{k}-u\right)+\frac{1}{h^{2}}\left(\frac{\partial{e}}{\partial{y}},\frac{\partial
e}{\partial
y}\right)+\frac{1}{h^{2}}\left(\frac{\partial{e}}{\partial{y}},\frac{\partial}{\partial
y}(w_{k}-u)\right)$
$\displaystyle-\left(\frac{1}{h_{k}^{2}}-\frac{1}{h^{2}}\right)\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial}{\partial
y}(w_{k}-u_{k})\right)-\left(\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial
e}{\partial
y}(w_{k}-u_{k})dy+\left(\frac{h^{\prime}}{h}-\frac{h_{k}^{\prime}}{h_{k}}\right)\int_{0}^{1}y\frac{\partial
u_{k}}{\partial y}(w_{k}-u_{k})dy\right)$
$\displaystyle+\frac{h^{\prime}}{h}e(\tau,1)(w_{k}(1)-u_{k}(1))+\left(\frac{h^{\prime}}{h}-\frac{h_{k}^{\prime}}{h_{k}}\right)u_{k}(\tau,1)(w_{k}(1)-u_{k}(1))$
$\displaystyle-\left(\text{Bi}\frac{b(\tau)}{m_{0}}\left(\frac{1}{h}-\frac{1}{h_{k}}\right)(w_{k}(0)-u_{k}(0))-\text{Bi}\;\text{H}\left(\frac{u(\tau,0)}{h}-\frac{u_{k}(\tau,0)}{h_{k}}\right)(w_{k}(0)-u_{k}(0))\right)=0.$

Therefore, we can write

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\lVert
e\rVert^{2}+\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle\leq\left\|\frac{\partial
e}{\partial\tau}\right\|\left\|u-w_{k}\right\|+\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|\left\|\frac{\partial}{\partial y}(u-w_{k})\right\|$
$\displaystyle+|h-h_{k}|\frac{h+h_{k}}{h^{2}h_{k}^{2}}\left\|\frac{\partial
u_{k}}{\partial y}\right\|\left\|\frac{\partial}{\partial
y}(w_{k}-u_{k})\right\|+\frac{h^{\prime}}{h}\left\|\frac{\partial e}{\partial
y}\right\|\left\|w_{k}-u_{k}\right\|$
$\displaystyle+\left|\frac{h^{\prime}}{h}-\frac{h^{\prime}_{k}}{h_{k}}\right|\left\|\frac{\partial
u_{k}}{\partial y}\right\|\left\|\frac{\partial}{\partial
y}(w_{k}-u_{k})\right\|+\frac{h^{\prime}}{h}|e(\tau,1)||(w_{k}(1)-u(1))+e(\tau,1)|$
$\displaystyle+\left|\frac{h^{\prime}}{h}-\frac{h^{\prime}_{k}}{h_{k}}\right||u_{k}(\tau,1)||(w_{k}(1)-u(1))+e(\tau,1)|$
$\displaystyle+\text{Bi}\frac{b^{*}}{m_{0}}\left|\frac{1}{h}-\frac{1}{h_{k}}\right||w_{k}(0)-u_{k}(0)|+\text{Bi}\;\text{H}\left|\frac{u(\tau,0)}{h}-\frac{u_{k}(\tau,0)}{h_{k}}\right||w_{k}(0)-u_{k}(0)|.$
(49)

To bound some terms on the right hand side in (49), we introduce the strictly
positive constant $c_{\ell}<\infty,\ell\in\\{1,2,\cdots,5\\}$. The value for
these constants is not explicitly written, but can be calculated. Before
proceeding further, we collect two useful estimates in Remark 5.1.

######  Remark 5.1.

There exist constants $c_{2},\;c_{5}>0$ such that

$\displaystyle(1)\;\;$
$\displaystyle\left|\frac{h^{\prime}}{h}-\frac{h^{\prime}_{k}}{h_{k}}\right|\leq
c_{2}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)$ $\displaystyle(2)\;\;$
$\displaystyle\left(\frac{u(0)}{h}-\frac{u_{k}(0)}{h_{k}}\right)=\frac{1}{h}(e(0))+\frac{u_{k}(0)}{h_{k}}(h_{k}-h)\leq
c_{5}(|e(0)|+|h-h_{k}|).$

Making use of Remark 5.1, (49) becomes

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\lVert
e\rVert^{2}+\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle\leq\left\|\frac{\partial
e}{\partial\tau}\right\|\left\|u-w_{k}\right\|+\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|\left\|\frac{\partial}{\partial y}(u-w_{k})\right\|$
$\displaystyle+c_{1}|h-h_{k}|\left\|\frac{\partial u_{k}}{\partial
y}\right\|\left(\left\|\frac{\partial}{\partial
y}(w_{k}-u)\right\|+\left\|\frac{\partial e}{\partial
y}\right\|\right)+\frac{h^{\prime}}{h}\left\|\frac{\partial e}{\partial
y}\right\|\left(\left\|w_{k}-u\right\|+\left\|e\right\|\right)$
$\displaystyle+c_{2}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\left\|\frac{\partial
u_{k}}{\partial y}\right\|\left(\left\|\frac{\partial}{\partial
y}(w_{k}-u)\right\|+\left\|\frac{\partial e}{\partial y}\right\|\right)$
$\displaystyle+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{h^{\prime}}{h}}|e(1)|(|w_{k}(1)-u(1)|+|e(1)|)$
$\displaystyle+c_{3}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)|u_{k}(\tau,1)|(|w_{k}(1)-u(1)|+|e(1)|)$
$\displaystyle+c_{4}\text{Bi}\frac{b^{*}}{m_{0}}\left|h-h_{k}\right|(|w_{k}(0)-u(0)|+|e(0)|)$
$\displaystyle+c_{5}\text{Bi}\;\text{H}(|e(0)|+|h-h_{k}|)(|w_{k}(0)-u(0)|+|e(0)|)=\sum_{\ell=1}^{9}I_{\ell}.$

We set $w_{k}:=I_{k}u$, where $I_{k}u$ is the Lagrange interpolation of $u$.
By using Lemma 3.1, Young’s inequality (8) and interpolation inequality (9),
we obtain the following estimates:

$\displaystyle I_{1}$ $\displaystyle:=\left\|\frac{\partial
e}{\partial\tau}\right\|\left\|u-w_{k}\right\|\leq\left\|\frac{\partial
e}{\partial\tau}\right\|\gamma_{1}k^{2}\|u\|_{H^{2}(0,1)}\leq\frac{1}{2}\left\|\frac{\partial
e}{\partial\tau}\right\|^{2}k^{2}+\frac{\gamma_{1}^{2}k^{2}}{2}\|u\|^{2}_{H^{2}(0,1)},$
$\displaystyle I_{2}$ $\displaystyle:=\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|\left\|\frac{\partial}{\partial
y}(u-w_{k})\right\|\leq\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|\gamma_{2}k\left\|u\right\|_{H^{2}(0,1)}\leq\frac{\xi}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+c_{\xi}\gamma_{2}^{2}k^{2}\frac{1}{h^{2}}\|u\|^{2}_{H^{2}(0,1)},$
$\displaystyle I_{3}$ $\displaystyle:=c_{1}|h-h_{k}|\left\|\frac{\partial
u_{k}}{\partial y}\right\|\left(\left\|\frac{\partial}{\partial
y}(w_{k}-u)\right\|+\left\|\frac{\partial e}{\partial y}\right\|\right)$
$\displaystyle\leq c_{1}|h-h_{k}|\left\|\frac{\partial u_{k}}{\partial
y}\right\|\left(\gamma_{2}k\left\|u\right\|_{H^{2}(0,1)}+\left\|\frac{\partial
e}{\partial y}\right\|\right)$
$\displaystyle\leq\rho|h-h_{k}|^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|^{2}}+c_{\rho}c_{1}^{2}\gamma_{2}^{2}k^{2}\left\|u\right\|^{2}_{H^{2}(0,1)}+c_{\hat{\rho}}c_{1}^{2}|h-h_{k}|^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|^{2}}h^{2}+\hat{\rho}\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|^{2}.$

We observe that if $u_{0}\in H^{1}(0,1)$, then it also holds that $u_{k}\in
H^{1}(S_{\hat{T}},V_{k})$. Hence, we can control terms like
$\left\|\frac{\partial u_{k}}{\partial y}\right\|^{2}$ via

$\displaystyle\max_{0\leq\tau\leq\hat{T}}\|u_{k}\|^{2}_{H^{1}(0,1)}+\int_{0}^{\hat{T}}\|u_{k}\|^{2}_{H^{2}(0,1)}d\tau\leq\hat{c}_{1}.$
(50)

Therefore, we get

$\displaystyle I_{3}$
$\displaystyle\leq\rho\hat{c}_{1}|h-h_{k}|^{2}+c_{\rho}c_{1}^{2}\gamma_{2}^{2}k^{2}\left\|u\right\|^{2}_{H^{2}(0,1)}+c_{\hat{\rho}}\hat{c}_{1}c_{1}^{2}|h-h_{k}|^{2}h^{2}+\hat{\rho}\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|^{2},$ $\displaystyle I_{4}$
$\displaystyle:=\frac{h^{\prime}}{h}\left\|\frac{\partial e}{\partial
y}\right\|\left(\left\|w_{k}-u\right\|+\left\|e\right\|\right)$
$\displaystyle\leq\frac{h^{\prime}}{h}\left\|\frac{\partial e}{\partial
y}\right\|\left(\gamma_{1}k^{2}\left\|u\right\|_{H^{2}(0,1)}+\left\|e\right\|\right)$
$\displaystyle\leq\zeta\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+c_{\zeta}(h^{\prime})^{2}\left(\gamma_{1}k^{2}\left\|u\right\|_{H^{2}(0,1)}+\left\|e\right\|\right)^{2}$
$\displaystyle\leq\zeta\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+2c_{\zeta}(h^{\prime})^{2}\left(\gamma_{1}^{2}k^{4}\left\|u\right\|^{2}_{H^{2}(0,1)}+\left\|e\right\|^{2}\right),$
$\displaystyle I_{5}$
$\displaystyle:=c_{2}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\left\|\frac{\partial
u_{k}}{\partial y}\right\|\left(\left\|\frac{\partial}{\partial
y}(w_{k}-u)\right\|+\left\|\frac{\partial e}{\partial y}\right\|\right)$
$\displaystyle\leq
c_{2}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|}\left(\gamma_{2}k\left\|u\right\|_{H^{2}(0,1)}+\left\|\frac{\partial
e}{\partial y}\right\|\right)$
$\displaystyle\leq\xi\left(|h-h_{k}|^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|^{2}}+|h^{\prime}-h^{\prime}_{k}|^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|^{2}}\right)+c_{\xi}c_{2}^{2}\gamma_{2}^{2}k^{2}\left\|u\right\|^{2}_{H^{2}(0,1)}$
$\displaystyle+\hat{\xi}\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+c_{\hat{\xi}}c_{2}^{2}h^{2}\left(|h-h_{k}|^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial
y}\right\|^{2}}+|h^{\prime}-h^{\prime}_{k}|^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\|\frac{\partial
u_{k}}{\partial y}\right\|^{2}}\right),$
$\displaystyle\leq\xi\hat{c}_{1}\left(|h-h_{k}|^{2}+|h^{\prime}-h^{\prime}_{k}|^{2}\right)+c_{\xi}c_{2}^{2}\gamma_{2}^{2}k^{2}\left\|u\right\|^{2}_{H^{2}(0,1)}$
$\displaystyle+\hat{\xi}\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+c_{\hat{\xi}}\hat{c}_{1}c_{2}^{2}h^{2}\left(|h-h_{k}|^{2}+|h^{\prime}-h^{\prime}_{k}|^{2}\right),$
$\displaystyle I_{6}$
$\displaystyle:=\frac{h^{\prime}}{h}|e(1)|(|w_{k}(1)-u(1)|+|e(1)|)$
$\displaystyle=\frac{h^{\prime}}{h}|e(1)|^{2}+\frac{h^{\prime}}{h}|e(1)||w_{k}(1)-u(1)|$
$\displaystyle=\frac{h^{\prime}}{h}|e(1)|^{2}+\frac{h^{\prime}}{h}\left(\frac{|e(1)|^{2}}{2}+\frac{|w_{k}(1)-u(1)|^{2}}{2}\right)$
$\displaystyle=\frac{3}{2}\frac{h^{\prime}}{h}|e(1)|^{2}+\frac{h^{\prime}}{h}\frac{|w_{k}(1)-u(1)|^{2}}{2}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{3}{2}\frac{h^{\prime}}{h}\hat{c}\left\|e\right\|^{2\theta}\left\|e\right\|_{H^{1}(0,1)}^{2(1-\theta)}+\frac{h^{\prime}}{2h}\left(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta}\right)^{2}\left\|u\right\|_{H^{2}(0,1)}^{2}}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{3}{2}\left(\frac{\xi}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+(\frac{\xi}{h^{2}}+c_{\xi}\hat{c}^{2}(h^{\prime})^{2})\left\|e\right\|^{2}\right)+\frac{h^{\prime}}{2h}\hat{c}^{2}\left(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta}\right)^{2}\left\|u\right\|_{H^{2}(0,1)}^{2},}$
$\displaystyle I_{7}$
$\displaystyle:=c_{3}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)|u_{k}(\tau,1)|(|w_{k}(1)-u(1)|+|e(1)|)$
$\displaystyle\leq
c_{3}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\hat{c}\left\|u_{k}\right\|^{1-\theta}\left\|\frac{\partial
u_{k}}{\partial y}\right\|^{\theta}(|w_{k}(1)-u(1)|+|e(1)|)$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}c_{3}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\hat{c}^{2}\left\|u_{k}\right\|_{H^{1}(0,1)}\left(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta}\right)\|u\|_{H^{2}(0,1)}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+c_{3}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\hat{c}^{2}\left\|u_{k}\right\|_{H^{1}(0,1)}\|e\|^{\theta}\left\|e\right\|_{H^{1}(0,1)}^{1-\theta}}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}c_{3}\hat{c}^{2}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\left\|u_{k}\right\|_{H^{1}(0,1)}\left(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta}\right)\|u\|_{H^{2}(0,1)}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+c_{3}\hat{c}^{2}(|h-h_{k}|+|h^{\prime}-h^{\prime}_{k}|)\left\|u_{k}\right\|_{H^{1}(0,1)}\|e\|^{\theta}\left\|e\right\|_{H^{1}(0,1)}^{1-\theta}}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\xi\left(|h-h_{k}|+|h^{\prime}-h_{k}^{\prime}|\right)^{2}\left\|u_{k}\right\|_{H^{1}(0,1)}^{2}+c_{\xi}\left(c_{3}\hat{c}^{2}(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta})\right)^{2}\|u\|_{H^{2}(0,1)}^{2}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\bar{\xi}\left(|h-h_{k}|+|h^{\prime}-h_{k}^{\prime}|\right)^{2}\left\|u_{k}\right\|_{H^{1}(0,1)}^{2}+c_{\bar{\xi}}c_{3}^{2}\hat{c}^{4}\|e\|^{2\theta}\left\|e\right\|_{H^{1}(0,1)}^{2(1-\theta)}}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2\xi\hat{c}_{1}\left(|h-h_{k}|^{2}+|h^{\prime}-h_{k}^{\prime}|^{2}\right)+c_{\xi}\left(c_{3}\hat{c}^{2}(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta})\right)^{2}\|u\|_{H^{2}(0,1)}^{2}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+2\bar{\xi}\hat{c}_{1}\left(|h-h_{k}|^{2}+|h^{\prime}-h_{k}^{\prime}|^{2}\right)+\hat{\xi}\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+\left(\frac{\hat{\xi}}{h^{2}}+c_{\hat{\xi}}c_{\bar{\xi}}^{2}c_{3}^{4}\hat{c}^{8}h^{2}\right)\|e\|^{2},}$
$\displaystyle I_{8}$
$\displaystyle:=c_{4}\text{Bi}\frac{b^{*}}{m_{0}}\left|h-h_{k}\right|(|w_{k}(0)-u(0)|+|e(0)|)$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}c_{4}\text{Bi}\frac{b^{*}}{m_{0}}\left|h-h_{k}\right|\left(\hat{c}(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta})\left\|u\right\|_{H^{2}(0,1)}+\hat{c}\left\|e\right\|^{\theta}\left\|e\right\|^{1-\theta}_{H^{1}(0,1)}\right)}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\xi|h-h_{k}|^{2}+c_{\xi}c_{4}^{2}\hat{c}^{2}\text{Bi}^{2}\frac{(b^{*})^{2}}{m_{0}^{2}}(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta})^{2}\|u\|_{H^{2}(0,1)}^{2}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\hat{\xi}|h-h_{k}|^{2}+c_{\hat{\xi}}c_{4}^{2}\hat{c}^{2}\text{Bi}^{2}\frac{(b^{*})^{2}}{m_{0}^{2}}\left\|e\right\|^{2\theta}\left\|e\right\|^{2(1-\theta)}_{H^{1}(0,1)}}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\xi|h-h_{k}|^{2}+c_{\xi}c_{4}^{2}\hat{c}^{2}\text{Bi}^{2}\frac{(b^{*})^{2}}{m_{0}^{2}}(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta})^{2}\|u\|_{H^{2}(0,1)}^{2}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\hat{\xi}|h-h_{k}|^{2}+\bar{\xi}\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+\left(\frac{\bar{\xi}}{h^{2}}+c_{\bar{\xi}}c_{\hat{\xi}}^{2}c_{4}^{4}\hat{c}^{4}\text{Bi}^{4}\frac{(b^{*})^{4}}{m_{0}^{4}}h^{2}\right)\left\|e\right\|^{2}.}$

By a similar calculation used to obtain the upper bounds on $I_{6}$ and
$I_{8}$, we get

$\displaystyle I_{9}$
$\displaystyle:=c_{5}\text{Bi}\;\text{H}(|e(0)|+|h-h_{k}|)(|w_{k}(0)-u(0)|+|e(0)|)$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{3}{2}\left(\frac{\xi}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+\left(\frac{\xi}{h^{2}}+c_{\xi}\hat{c}^{2}c_{5}^{2}\text{Bi}^{2}\text{H}^{2}h^{2}\right)\left\|e\right\|^{2}\right)+c_{\xi}c_{5}^{2}\text{Bi}^{2}\text{H}^{2}\hat{c}^{2}(\gamma_{1}k^{2}+\gamma_{3}k^{1+\theta})^{2}\left\|u\right\|_{H^{2}(0,1)}^{2}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+(\xi+\hat{\xi})|h-h_{k}|^{2}+\bar{\xi}\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+\left(\frac{\bar{\xi}}{h^{2}}+c_{\bar{\xi}}c_{\hat{\xi}}^{2}c_{5}^{4}\hat{c}^{4}\text{Bi}^{4}\text{H}^{4}h^{2}\right)\left\|e\right\|^{2}.}$

Finally, we are led to the following structural inequality:

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\lVert
e\rVert^{2}+A_{1}\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle\leq
A_{2}k^{2}+A_{3}\left\|e\right\|^{2}+A_{4}|h-h_{k}|^{2}+A_{5}|h^{\prime}-h_{k}^{\prime}|^{2},$
(51)

where

$\displaystyle A_{1}$
$\displaystyle:=\frac{1}{L^{2}}\left(1-\frac{5}{2}\xi-\hat{\rho}-\zeta-2\hat{\xi}-2\bar{\xi}\right),$
$\displaystyle A_{2}$
$\displaystyle:=\|u\|_{H^{2}(0,1)}\left(\frac{\gamma_{1}^{2}}{2}+\frac{1}{h_{0}^{2}}c_{\xi}\gamma_{2}^{2}+c_{\rho}c_{1}^{2}\gamma_{2}^{2}+2c_{\zeta}\|h^{\prime}\|_{\infty}^{2}\gamma_{1}^{2}+c_{\xi}c_{2}^{2}\gamma_{2}^{2}+\frac{\|h^{\prime}\|_{\infty}}{2h_{0}}\hat{c}^{2}(\gamma_{1}+\gamma_{3})^{2}\right.$
$\displaystyle\left.+c_{\xi}c_{3}^{2}\hat{c}^{4}(\gamma_{1}+\gamma_{3})^{2}+c_{\xi}c_{4}^{2}\hat{c}^{2}\text{Bi}^{2}\frac{(b^{*})^{2}}{m_{0}^{2}}(\gamma_{1}+\gamma_{3})^{2}+c_{\xi}c_{5}^{2}\hat{c}^{2}\text{Bi}^{2}\;\text{H}^{2}(\gamma_{1}+\gamma_{3})^{2}\right),$
$\displaystyle A_{3}$
$\displaystyle:=2c_{\zeta}\|h^{\prime}\|_{\infty}^{2}+\frac{1}{h_{0}^{2}}\left(3\xi+\hat{\xi}+2\bar{\xi}\right)+\frac{3}{2}c_{\xi}\hat{c}^{2}\|h^{\prime}\|^{2}_{\infty}+c_{\hat{\xi}}c_{\bar{\xi}}^{2}c_{3}^{4}\hat{c}^{8}\|h\|^{2}_{\infty}+c_{\bar{\xi}}c_{\hat{\xi}}^{2}\hat{c}^{4}\text{Bi}^{4}\left(c_{4}^{4}\frac{(b^{*})^{4}}{m_{0}^{4}}+c_{5}^{4}\rm{H}^{4}\right)\|h\|^{2}_{\infty}$
$\displaystyle+c_{\xi}\hat{c}^{2}c_{5}^{2}\text{Bi}^{2}\;\text{H}^{2}\|h\|^{2}_{\infty},$
$\displaystyle A_{4}$
$\displaystyle:=2(\xi+\hat{\xi})+\hat{c}_{1}(\rho+\xi)+2\hat{c}_{1}(\bar{\xi}+\xi)+c_{\hat{\rho}}\hat{c}_{1}c_{1}^{2}\|h\|_{\infty}^{2}+c_{2}^{2}\hat{c}_{1}c_{\hat{\xi}}\|h\|_{\infty}^{2},$
$\displaystyle A_{5}$
$\displaystyle:=3\xi\hat{c}_{1}+2\bar{\xi}\hat{c}_{1}+c_{2}^{2}\hat{c}_{1}c_{\hat{\xi}}\|h\|_{\infty}^{2}.$

From (19) and (23), we get for all $\tau\in(0,\hat{T})$ the inequality

$\displaystyle|h^{\prime}(\tau)-h^{\prime}_{k}(\tau)|$ $\displaystyle\leq
A_{0}|e(1)|+\frac{1}{m_{0}}|\sigma(h(\tau))-\sigma(h_{k}(\tau))|$
$\displaystyle\leq
A_{0}\hat{c}\left(\eta\left\|e(\tau)\right\|_{H^{1}(0,1)}+c_{\eta}\|e(\tau)\|\right)+\frac{\mathcal{L}}{m_{0}}|h(\tau)-h_{k}(\tau)|.$

Thus, this leads to

$\displaystyle|h^{\prime}-h^{\prime}_{k}|^{2}\leq
3\left(A_{0}^{2}\hat{c}^{2}\eta^{2}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+A_{0}^{2}\hat{c}^{2}(\eta^{2}+c_{\eta}^{2})\|e\|^{2}+\frac{\mathcal{L}^{2}}{m_{0}^{2}}|h-h_{k}|^{2}\right).$
(52)

Using (52) in (51), we infer that

$\displaystyle\frac{d}{d\tau}\lVert
e\rVert^{2}+(A_{1}-3A_{0}^{2}\hat{c}^{2}\eta^{2}A_{5})\left\|\frac{\partial
e}{\partial y}\right\|^{2}\leq A_{2}k^{2}$
$\displaystyle+A_{6}\left\|e\right\|^{2}+\left(A_{4}+3A_{5}\frac{\mathcal{L}^{2}}{m_{0}^{2}}\right)|h-h_{k}|^{2},$
(53)

where $A_{6}:=A_{3}+3A_{0}^{2}\hat{c}^{2}(\eta^{2}+c_{\eta}^{2})A_{5}$. We
choose $\xi>0,\;\hat{\rho}>0,\;\bar{\xi}>0,\;\zeta>0$, $\hat{\xi}>0$, and
$\eta>0$ sufficiently small such that
$\zeta_{1}:=A_{1}-3A_{0}^{2}\hat{c}^{2}\eta^{2}A_{5}>0$. Applying Gronwall’s
inequality (see e.g. Appendix B in [9]) gives the following upper bounds:

$\displaystyle\|e(\tau)\|^{2}$ $\displaystyle\leq
e^{\displaystyle\int_{0}^{\tau}A_{6}ds}\left(\|e(0)\|^{2}+\int_{0}^{\tau}\left(A_{2}k^{2}+\left(A_{4}+3A_{5}\frac{\mathcal{L}^{2}}{m_{0}^{2}}\right)|h(s)-h_{k}(s)|^{2}\right)ds\right)$
$\displaystyle\leq
c_{6}(A_{0},A_{3},A_{5},\hat{T})\left(k^{4}\|u_{0}\|^{2}_{H^{2}(0,1)}+A_{2}k^{2}\tau+\left(A_{4}+3A_{5}\frac{\mathcal{L}^{2}}{m_{0}^{2}}\right)\int_{0}^{\tau}|h(s)-h_{k}(s)|^{2}ds\right)$
$\displaystyle\leq
c_{6}(A_{0},A_{3},A_{4},A_{5},\mathcal{L},\hat{T})\left(k^{4}+k^{2}\hat{T}+\|h-h_{k}\|^{2}_{L^{2}(S_{\hat{T}})}\right).$

Thus, we obtain

$\displaystyle\displaystyle\max_{0\leq\tau\leq\hat{T}}\|e(\tau)\|^{2}\leq
c_{6}\left(k^{2}+\|h-h_{k}\|^{2}_{L^{2}(S_{\hat{T}})}\right).$

By using Young’s inequality together with (52), we get the following
relations:

$\displaystyle\frac{d}{d\tau}(|h-h_{k}|^{2})$
$\displaystyle=2(h-h_{k})(h^{\prime}-h_{k}^{\prime})$
$\displaystyle\leq|h-h_{k}|^{2}+|h^{\prime}-h_{k}^{\prime}|^{2}$
$\displaystyle\leq
C|h-h_{k}|^{2}+3A_{0}^{2}\hat{c}^{2}\eta^{2}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+3A_{0}^{2}\hat{c}^{2}(\eta^{2}+c_{\eta}^{2})\|e\|^{2},$ (54)

where $C:=1+3\mathcal{L}^{2}/m_{0}^{2}.$
Let $\delta>0$ be any positive real number. Adding
$\delta\frac{d}{d\tau}|h-h_{k}|^{2}$ on both sides of (53) and using (54)
yields

$\displaystyle\frac{d}{d\tau}\left(\lVert
e\rVert^{2}+\delta|h-h_{k}|^{2}\right)+(\zeta_{1}-3\delta\hat{c}^{2}A_{0}^{2}\eta^{2})\left\|\frac{\partial
e}{\partial y}\right\|^{2}\leq A_{2}k^{2}$ $\displaystyle+(A_{6}+3\delta
A_{0}^{2}\hat{c}^{2}(\eta^{2}+c_{\eta}^{2}))\|e\|^{2}$
$\displaystyle+\left(A_{4}+3A_{5}\frac{\mathcal{L}^{2}}{m_{0}^{2}}+\delta
C\right)|h-h_{k}|^{2}.$ (55)

We choose $\eta>0$ in such a way that
$(\zeta_{1}-3\delta\hat{c}^{2}A_{0}^{2}\eta^{2})>0$. Then it exists a constant
$A_{7}>0$ such that

$\displaystyle\frac{d}{d\tau}\left(\lVert
e\rVert^{2}+\delta|h-h_{k}|^{2}\right)\leq
A_{2}k^{2}+A_{7}(\|e\|^{2}+\delta|h-h_{k}|^{2}).$ (56)

Gronwall’s inequality applied to (56) for the quantity
$\|e\|^{2}+\delta|h-h_{k}|^{2}$ gives the estimate

$\displaystyle\|e\|^{2}+\delta|h-h_{k}|^{2}\leq ck^{2}.$ (57)

Integrating (53) from $0$ to $\hat{T}$ and using (57) yields

$\displaystyle\int_{0}^{\hat{T}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}d\tau\leq c_{7}k^{2}.$ (58)

Integrating (52) from $0$ to $\hat{T}$ and using (57) and (58) gives the
estimate

$\displaystyle\|h^{\prime}-h_{k}^{\prime}\|^{2}\leq ck^{2},$

which completes the proof of Theorem 5.1. ∎

######  Theorem 5.2.

(A posteriori error estimate) Assume (A1)–(A5) hold. Additionally, take
$u_{0}\in H^{2}(0,1)$. Let $(u,h)$ and $(u_{k},h_{k})$ be the corresponding
weak solutions to the problem $(P)$ and $(P_{d})$ in the sense of Definition
4.1 and Definition 4.2, respectively. Then there exist
$0<\tilde{T}\leq\hat{T}$ and positive constants $c_{1},\;c_{2},\,c_{3}$
(independent of $k$ and $u$) such that for all $\tau\in
S_{\tilde{T}}:=(0,\tilde{T})$ the following inequality holds:

$\displaystyle\|u-u_{k}\|_{L^{2}(0,1)}+c_{1}|h-h_{k}|^{2}+c_{2}\int_{0}^{\tau}\left\|\frac{\partial}{\partial
x}(u-u_{k})\right\|^{2}ds$ $\displaystyle\leq
c_{3}\left(|h(0)-h_{k}(0)|^{2}+\sum_{i=0}^{N-2}k_{i}^{2}\left\\{\|R(u_{k})\|_{L^{2}(S_{\tilde{T}},L^{2}(I_{i}))}^{2}+k_{i}^{2}\|u_{0}\|_{H^{2}(I_{i})}^{2}\right\\}\right),$
(59)

where the residual $R(u_{k})$ is defined by

$\displaystyle R(u_{k}):=\frac{h^{\prime}_{k}}{h_{k}}y\frac{\partial
u_{k}}{\partial y}+\frac{1}{h_{k}}{\rm Bi}\left(\frac{b(\tau)}{m_{0}}-{\rm
H}u_{k}(\tau,0)\right)-\frac{h^{\prime}_{k}}{h_{k}}u_{k}(\tau,1)-\frac{\partial
u_{k}}{\partial\tau}.$ (60)

###### Proof.

Let $e:=u-u_{k}$ be the pointwise error. Using the weak formulation (18), we
can write

$\displaystyle\left(\frac{\partial
e}{\partial\tau},v\right)+\frac{1}{h^{2}}\left(\frac{\partial e}{\partial
y},\frac{\partial v}{\partial y}\right)=$
$\displaystyle\left[\left(\frac{\partial
u}{\partial\tau},v\right)+\frac{1}{h^{2}}\left(\frac{\partial u}{\partial
y},\frac{\partial v}{\partial y}\right)\right]-\left[\left(\frac{\partial
u_{k}}{\partial\tau},v\right)+\frac{1}{h^{2}}\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v}{\partial y}\right)\right]$
$\displaystyle=$ $\displaystyle\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial
u}{\partial
y}vdy+\frac{1}{h}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u(\tau,0)\right)v(0)-\frac{h^{\prime}}{h}u(\tau,1)v(1)$
$\displaystyle-\left[\left(\frac{\partial
u_{k}}{\partial\tau},v\right)+\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v}{\partial
y}\right)+\left(\frac{1}{h^{2}}-\frac{1}{h_{k}^{2}}\right)\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v}{\partial y}\right)\right]$ (61)

for all $v\in H^{1}(0,1)$. Inserting (60) into (61) yields

$\displaystyle\left(\frac{\partial
e}{\partial\tau},v\right)+\frac{1}{h^{2}}\left(\frac{\partial e}{\partial
y},\frac{\partial v}{\partial y}\right)=$
$\displaystyle\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial u}{\partial
y}vdy+\frac{1}{h}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u(\tau,0)\right)v(0)-\frac{h^{\prime}}{h}u(\tau,1)v(1)$
$\displaystyle-\left(\frac{1}{h^{2}}-\frac{1}{h_{k}^{2}}\right)\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v}{\partial
y}\right)-\frac{h^{\prime}_{k}}{h_{k}}\int_{0}^{1}y\frac{\partial
u_{k}}{\partial
y}vdy-\frac{1}{h_{k}}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u_{k}(\tau,0)\right)v(0)$
$\displaystyle+\frac{h^{\prime}_{k}}{h_{k}}u_{k}(\tau,1)v(1)+\left[\int_{0}^{1}R(u_{k})vdy-\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v}{\partial y}\right)\right],$ (62)

where $R(u_{k})$ is the residual quantity defined in (60). Since $u_{k}\in
V_{k}$, we have that $\displaystyle\frac{\partial^{2}u_{k}}{\partial y^{2}}=0$
on each $I_{i}:=(y_{i},y_{i+1})$. The term

$\displaystyle\int_{0}^{1}R(u_{k})vdy-\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial v}{\partial y}\right)$

becomes after integration by part

$\displaystyle\sum_{i=0}^{N-2}\left\\{\int_{y_{i}}^{y_{i+1}}R(u_{k})vdy-\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y}(y_{i+1})v(y_{i+1})-\frac{\partial u_{k}}{\partial
y}(y_{i})v(y_{i})\right)\right\\}.$

We also get from (22)

$\displaystyle\sum_{i=0}^{N-2}\left\\{\int_{y_{i}}^{y_{i+1}}R(u_{k})v_{k}dy-\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y}(y_{i+1})v_{k}(y_{i+1})-\frac{\partial u_{k}}{\partial
y}(y_{i})v_{k}(y_{i})\right)\right\\}=0$ (63)

for all $v_{k}\in V_{k}$. Adding (63) to (62) while taking $v=e\in H^{1}(0,1)$
and $v_{k}=I_{k}e\in V_{k}$ gives

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\lVert
e\rVert^{2}+\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle=\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial u}{\partial
y}edy+\frac{1}{h}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u(\tau,0)\right)e(0)$
$\displaystyle-\frac{h^{\prime}}{h}u(\tau,1)e(1)-\left(\frac{1}{h^{2}}-\frac{1}{h_{k}^{2}}\right)\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial e}{\partial
y}\right)-\frac{h^{\prime}_{k}}{h_{k}}\int_{0}^{1}y\frac{\partial
u_{k}}{\partial y}edy$
$\displaystyle-\frac{1}{h_{k}}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u_{k}(\tau,0)\right)e(0)+\frac{h^{\prime}_{k}}{h_{k}}u_{k}(\tau,1)e(1)$
$\displaystyle+\sum_{i=0}^{N-2}\Biggl{\\{}\int_{y_{i}}^{y_{i+1}}R(u_{k})(e-I_{k}e)dy$
$\displaystyle-\frac{1}{h_{k}^{2}}\Biggl{(}\frac{\partial u_{k}}{\partial
y}(y_{i+1})(e-I_{k}e)(y_{i+1})-\frac{\partial u_{k}}{\partial
y}(y_{i})(e-I_{k}e)(y_{i})\Biggr{)}\Biggr{\\}}$
$\displaystyle=\sum_{i=1}^{5}I_{i},$

where

$\displaystyle I_{1}$
$\displaystyle:=\frac{h^{\prime}}{h}\int_{0}^{1}y\frac{\partial u}{\partial
y}edy-\frac{h^{\prime}_{k}}{h_{k}}\int_{0}^{1}y\frac{\partial u_{k}}{\partial
y}edy,$ $\displaystyle I_{2}$
$\displaystyle:=\frac{1}{h}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u(\tau,0)\right)e(0)-\frac{1}{h_{k}}\text{Bi}\left(\frac{b(\tau)}{m_{0}}-\text{H}u_{k}(\tau,0)\right)e(0),$
$\displaystyle I_{3}$
$\displaystyle:=\frac{h^{\prime}_{k}}{h_{k}}u_{k}(\tau,1)e(1)-\frac{h^{\prime}}{h}u(\tau,1)e(1),$
$\displaystyle I_{4}$
$\displaystyle:=-\left(\frac{1}{h^{2}}-\frac{1}{h_{k}^{2}}\right)\left(\frac{\partial
u_{k}}{\partial y},\frac{\partial e}{\partial y}\right),$ $\displaystyle
I_{5}$
$\displaystyle:=\sum_{i=0}^{N-2}\left\\{\int_{y_{i}}^{y_{i+1}}R(u_{k})(e-I_{k}e)dy-\frac{1}{h_{k}^{2}}\left(\frac{\partial
u_{k}}{\partial y}(y_{i+1})(e-I_{k}e)(y_{i+1})-\frac{\partial u_{k}}{\partial
y}(y_{i})(e-I_{k}e)(y_{i})\right)\right\\}.$

By using (26) together with Cauchy-Schwarz and Young’s inequality, we obtain

$\displaystyle|I_{1}|$
$\displaystyle\leq\frac{h^{\prime}}{h}\left\|\frac{\partial e}{\partial
y}\right\|\left\|e\right\|+\left|\frac{h^{\prime}}{h}-\frac{h^{\prime}_{k}}{h_{k}}\right|\|e\|$
$\displaystyle\leq\left(\frac{\xi}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+c_{\xi}\|h^{\prime}\|_{\infty}^{2}\left\|e\right\|^{2}\right)+\xi\|e\|^{2}+2c_{\xi}\left(|h-h_{k}|^{2}+|h^{\prime}-h_{k}^{\prime}|^{2}\right).$
(64) $\displaystyle|I_{2}|$
$\displaystyle\leq\text{Bi}\;\frac{b(\tau)}{m_{0}}\frac{1}{hh_{k}}|h-h_{k}||e(0)|+\text{Bi}\;\text{H}\left|\frac{u(\tau,0)}{h}-\frac{u_{k}(\tau,0)}{h_{k}}\right||e(0)|$
$\displaystyle\leq\left(\text{Bi}\frac{b^{*}}{m_{0}}\frac{1}{L^{2}}\hat{c}+\text{Bi}\;\text{H}\hat{c}\right)|h-h_{k}|\|e\|^{1-\theta}\left\|\frac{\partial
e}{\partial
y}\right\|^{\theta}+c_{2}\text{Bi}\;\text{H}\hat{c}\|e\|^{2(1-\theta)}\left\|\frac{\partial
e}{\partial y}\right\|^{2\theta}$ $\displaystyle\leq\bar{\xi}|h-h_{k}|^{2}+\xi
c_{\bar{\xi}}\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+\tilde{c}^{4}c_{\xi}c_{\bar{\xi}}h^{2}\|e\|^{2}+\frac{\xi}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|^{2}+\tilde{c}_{1}^{2}c_{\xi}h^{2}\|e\|^{2},$ (65)

where

$\theta=\frac{1}{2},\;\;\tilde{c}:=\left(\text{Bi}\frac{b^{*}}{m_{0}}\frac{1}{L^{2}}\hat{c}+\text{Bi}\;\text{H}\hat{c}\right)\;\;\;\;\text{and}\;\;\;\tilde{c}_{1}:=c_{2}\text{Bi}\;\text{H}\hat{c}.$
$\displaystyle|I_{3}|$
$\displaystyle\leq\left|\frac{h^{\prime}}{h}-\frac{h^{\prime}_{k}}{h_{k}}\right||e(1)|+\frac{h_{k}^{\prime}}{h_{k}}|e(1)|^{2}$
$\displaystyle\leq
2\bar{\xi}\left(|h-h_{k}|^{2}+|h^{\prime}-h_{k}^{\prime}|^{2}\right)+\xi\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial
y}\right\|^{2}+c_{\bar{\xi}}c_{\xi}c_{3}^{4}\hat{c}^{4}\|e\|^{2}h^{2}+c\left(\frac{\xi}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|^{2}+c_{\xi}\|e\|^{2}\right).$ (66)
$\displaystyle|I_{4}|$
$\displaystyle\leq|h-h_{k}|\frac{h+h_{k}}{h^{2}h_{k}^{2}}\left\|\frac{\partial
u_{k}}{\partial y}\right\|\left\|\frac{\partial e}{\partial y}\right\|$
$\displaystyle\leq\xi|h-h_{k}|^{2}+c_{\xi}c^{2}(h_{0},L)\frac{1}{h^{2}}\left\|\frac{\partial
e}{\partial y}\right\|^{2}.$ (67)

To bound $|I_{5}|$ from above, we use the fact that $I_{k}e$ is the Lagrange
interpolant of $e$ with the property
$(e-I_{k}e)(y_{i})=0,\;i\in\\{0,1,2,\cdots,N-1\\}.$ We have

$\displaystyle|I_{5}|$
$\displaystyle\leq\sum_{i=0}^{N-2}\int_{y_{i}}^{y_{i+1}}R(u_{k})(e-I_{k}e)dy$
$\displaystyle\leq\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}\|e-I_{k}e\|_{L^{2}(I_{i})}$
$\displaystyle\leq\tilde{c}\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}k_{i}\left\|\frac{\partial
e}{\partial y}\right\|_{L^{2}(I_{i})}$
$\displaystyle\leq\tilde{c}\left(\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}\right)^{\frac{1}{2}}\left(\sum_{i=0}^{N-2}\left\|\frac{\partial
e}{\partial y}\right\|_{L^{2}(I_{i})}^{2}\right)^{\frac{1}{2}}$
$\displaystyle=\tilde{c}\left(\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}\right)^{\frac{1}{2}}\left\|\frac{\partial
e}{\partial y}\right\|_{L^{2}(0,1)}.$

By using Young’s inequality, we obtain

$\displaystyle|I_{5}|\leq\frac{\xi}{h^{2}}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+c_{\xi}\tilde{c}^{2}h^{2}\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}.$
(68)

It follows from (64)–(68) that for all $\xi,\;\bar{\xi}>0$, there exist
positive constants $K_{1},\;K_{2},\;K_{3}$ and $K_{4}$ such that

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\lVert
e\rVert^{2}+\frac{1}{h^{2}}\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle\leq K_{1}\|e\|^{2}+K_{2}|h-h_{k}|^{2}$
$\displaystyle+\frac{1}{h^{2}}K_{3}\left\|\frac{\partial e}{\partial
y}\right\|^{2}+K_{4}\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}.$

Let $\delta>0$ be a fixed, sufficiently small. Adding
$\frac{\delta}{2}\frac{d}{d\tau}|h-h_{k}|^{2}$ on both sides and using (54)
yields

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\left(\lVert
e\rVert^{2}+\delta|h-h_{k}|^{2}\right)+\frac{1}{L^{2}}(1-K_{3}-3\delta
A_{0}^{2}\eta)\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle\leq K_{1}\|e\|^{2}+K_{2}|h-h_{k}|^{2}+3\delta
A_{0}^{2}c_{\eta}\|e\|^{2}$
$\displaystyle+C\delta\|h-h_{k}\|^{2}+K_{4}\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}.$

We choose $\xi>0,\;\bar{\xi}>0$ and $\eta>0$ in such a way that
$1-K_{3}-3\delta A_{0}^{2}\eta\geq 0$. Then it exists $K_{5}>0$ such that

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\left(\lVert
e\rVert^{2}+\delta|h-h_{k}|^{2}\right)+\frac{1}{L^{2}}(1-K_{3}-3\delta
A_{0}^{2}\eta)\left\|\frac{\partial e}{\partial y}\right\|^{2}$
$\displaystyle\leq K_{5}(\|e\|^{2}+\delta|h-h_{k}|^{2})+$
$\displaystyle+K_{4}\sum_{i=0}^{N-2}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}.$
(69)

Applying Gronwall’s inequality to (69) for the quantity $\lVert
e\rVert^{2}+\delta|h-h_{k}|^{2}$ and using the initial condition

$\|e(0)\|_{L^{2}(0,1)}^{2}=\sum_{i=0}^{N-2}\|e(0)\|_{L^{2}(I_{i})}^{2}\leq
k_{i}^{4}\|u_{0}\|_{H^{2}(I_{i})}^{2},$

it exists a constant $c(\tilde{T},L)$ such that

$\displaystyle\lVert e\rVert^{2}+\delta|h(\tau)-h_{k}(\tau)|^{2}\leq
c(\tilde{T},L)\left(|h(0)-h_{k}(0)|^{2}+k_{i}^{4}\|u_{0}\|_{H^{2}(I_{i})}^{2}+\sum_{i=0}^{N-2}\int_{0}^{\tau}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}ds\right).$
(70)

By integrating (69) on $(0,\tau)$ and by using (70), it exists another
constant $c(\tilde{T},L)>0$ such that the following inequality holds:

$\displaystyle\int_{0}^{\tau}\left\|\frac{\partial}{\partial
x}(u-u_{k})\right\|^{2}ds\leq
c(\tilde{T},L)\left(|h(0)-h_{k}(0)|^{2}+k_{i}^{4}\|u_{0}\|_{H^{2}(I_{i})}^{2}+\sum_{i=0}^{N-2}\int_{0}^{\tau}\|R(u_{k})\|_{L^{2}(I_{i})}^{2}k_{i}^{2}ds\right).$

This concludes the proof of Theorem 5.2. ∎

##  6 Numerical illustrations

In this section, we firstly present our simulation results for both the dense
and foam rubber. The difference in the two cases is incorporated in the choice
of parameters. To approximate numerically the weak solution to (22)–(25), we
use the method of lines; for more details see, for instance, [20]. Firstly,
the model equations are discretized in space by means of the finite element
method. The resulting time-dependent system of ordinary differential equations
is tackled via the solver odeint in Python; see [21] for details on Python and
[13] for details on the solver. We refer the reader to see our previous work
[26] for the laboratory experiments, numerical method and simulation results
where we investigated the parameter space by exploring eventual effects of the
choice of parameters on the overall diffusants penetration process.
We take as observation time $T_{f}=40$ minutes for the final time with time
step $\Delta t=1/1000$ minutes. We choose the number of space discretization
points $N$ to be $100$. The values of parameters are taken to be $s_{0}=0.01$
(mm), $m_{0}$ = 0.1 (gram/mm3) and $b=1$ (gram/mm3). We take the value
$3.66\times 10^{-4}$ (mm2/min) for the diffusion constant $D$ [23], $0.564$
(mm/min) for absorption rate $\beta$ [31] and $2.5$ for Henry’s constant H
[5]. For the dense rubber, we choose $\sigma(s(t))=s(t)/10$ (gram/mm3) and
$a_{0}=500$ (mm4/sec/gram) while we choose $\sigma(s(t))=s(t)/50$ (gram/mm3)
and $a_{0}=2000$ (mm4/sec/gram) for the foam rubber case.

Figure 2: Dense rubber case. Left: Concentration profile of diffusant. Right:
Position of the moving boundary.

Figure 3: Foam rubber case. Left: Concentration profile of diffusants. Right:
Position of the moving boundary.

In Figure 2 and Figure 3 we show the concentration profile of the penetrating
diffusant, and respectively, the position of the moving boundary for the dense
rubber and foam rubber respectively. Comparing the diffusant concentration
profile in Figure 2 and Figure 3, we notice in both cases that, within a short
time of release of diffusant from its initial position, the diffusant quickly
enters the rubber from the left boundary and then starts diffusing inside
displacing a penetration front. In bothe Figure 2 and Figure 3, we compare
simulation results against experimental data for the position of moving
boundary. Both plots show a good agreement between model and experiment.

Finally, we wish to point out that the order of convergence of our FEM scheme
is consistent with the estimates stated in (46). As we are not aware of an
exact solution to (22)–(25), we compute the finite element approximation of
our weak solution on a fine mesh (say, with $640$ nodes) and denote it by
$u_{\tilde{k}}$. We use this $u_{\tilde{k}}$ as the reference solution for
computing the errors and convergence orders. We make use of the discrete
$\ell^{2}(Q(\hat{T}))$ norm which we denote here as

$\displaystyle
e(k_{i}):=\left\|u_{\tilde{k}}(\tau,y)-u_{k_{i}}(\tau,y)\right\|_{L^{2}(S_{\hat{T}},L^{2}(0,1))}=\left(\Delta\tau
k_{i}\sum_{j=0}^{N_{t}}\sum_{\ell=0}^{N-1}|u_{\tilde{k}}(\tau_{j},y_{\ell})-u_{k_{i}}(\tau_{j},y_{\ell})|^{2}\right)^{\frac{1}{2}}.$

Here $\Delta\tau$ is the uniform size of the $N_{t}+1$ time steps, while
$\\{k_{1},k_{2},k_{3},\cdots\\}$ with $k_{i}>k_{i+1}$ for
$i\in\\{1,2,\cdots\\}$ is a finite collection of the different mesh sizes used
in the computations.

We determine the convergence order based on any two consecutive calculations
of discrete errors using two different mesh sizes. To this end, we perform the
computations on a sequence of grids with mesh size $k$ that are halved in each
step. Thus, we use the following formula to compute the convergence order $r$:

$\displaystyle r:=\log_{2}\left(\frac{e(k_{i})}{e(k_{i+1})}\right).$

Figure 4:  Convergence order when time step size $\Delta t=10^{-4}$ is fixed.
Dash lines are lines of slope 1. Left: Log log scale plot of error on the
boundary $\|h-h_{k}\|_{L^{2}(S_{\hat{T}})}$ (circles) and
$\|h^{\prime}-h_{k}^{\prime}\|_{L^{2}(S_{\hat{T}})}$ (diamonds). Right: Log
log scale plot of error on the concentration
$\|u-u_{k}\|_{L^{2}(S_{\hat{T}},L^{2}(0,1))}$ (triangles) and
$\|u-u_{k}\|_{L^{2}(S_{\hat{T}},H^{1}(0,1))}$ (squares).

We show in Figure 4 the computed convergence order for the approximation of
the moving boundary position and of the concentration profile. This is done in
various norms for $N=20,40,80,160$, and $320$. These numerical results are in
agreement with the convergence order proven in Section 5.

##  7 Conclusion

The goal of this work was to analyze the errors produced by a semi-discrete
finite element approximation of the weak solution of moving boundary problem
modeling the penetration of diffusants into rubber. We obtained the a priori
error estimate (46) for the diffusant concentration profile as well as for the
position and speed of the moving boundary. The convergence rate is of order of
$\mathcal{O}(1)$ – the deviation from optimality is due to the nonlinear
coupling produced by the presence of the unknown moving boundary.
Additionally, we obtained the a posteriori error (59). Finally, we illustrated
numerically the basic output of our model. It turns out that results are in
the expected experimental range and they can be obtained in practice using
convergence rates closed to the theoretical ones.

## Acknowledgements

The authors acknowledge fruitful discussions with U. Giese, N. Kröger, R.
Meyer (Deutsches Institut für Kautschuktechnologie, Hannover, Germany), T.
Aiki (Japan Women’s University, Tokyo, Japan), and K. Kumazaki (Nagasaki
University, Japan) about the modeling, mathematical analysis, and simulation
of rubber-based materials exposed to environmental conditions. The work of
S.N. and A.M. is financed partly by the Swedish Research Council’s project
”Homogenization and dimension reduction of thin heterogeneous layers”, grant
nr. VR 2018-03648. A.M. also thanks the Knowledge Foundation for the grant KK
2019-0213, which led to the formulation of this problem setting.

## References

  * [1] R. A. Adams and J. F. Fournier.  Sobolev Spaces.  Elsevier, 2003.
  * [2] T. Aiki, N. H. Kröger, and A. Muntean.  A macro-micro elasticity-diffusion system modeling absorption-induced swelling in rubber foams–proof of the strong solvability.  Quarterly of Applied Mathematics, 79:545–579, 2021.
  * [3] T. Aiki and A. Muntean.  A free-boundary problem for concrete carbonation: Front nucleation and rigorous justification of the $\sqrt{t}$-law of propagation.  Interfaces and Free Boundaries, 15:167–180, 2013.
  * [4] V. Alexiades.  Mathematical Modeling of Melting and Freezing Processes.  CRC Press, 1992.
  * [5] M. Böhm, J. Devinny, F. Jahani, and G. Rosen.  On a moving-boundary system modeling corrosion in sewer pipes.  Applied Mathematics and Computation, 92(2-3):247–269, 1998.
  * [6] J.H. Brandts, S. Korotov, and M. Krizek.  The discrete maximum principle for linear simplicial finite element approximations of a reaction–diffusion problem.  Linear Algebra and its Applications, 99:2344–2357, 2008.
  * [7] A. Caboussat and J. Rappaz.  Analysis of a one-dimensional free boundary flow problem.  Numerische Mathematik, 101(1):67–86, 2005.
  * [8] C. Chainais-Hillairet, B. Merlet, and A. Zurek.  Convergence of a finite volume scheme for a parabolic system with a free boundary modeling concrete carbonation.  ESAIM: Mathematical Modelling and Numerical Analysis, 52(2):457–480, 2018.
  * [9] L. C. Evans.  Partial Differential Equations.  Graduate studies in mathematics. American Mathematical Society, 2010.
  * [10] S. C. Gupta.  The Classical Stefan Problem: Basic Concepts, Modelling and Analysis with Quasi-analytical Solutions and Methods, volume 45.  Elsevier, 2017.
  * [11] J. G. Heywood and R. Rannacher.  Finite element approximation of the nonstationary Navier–Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization.  SIAM Journal on Numerical Analysis, 19(2):275–311, 1982.
  * [12] E. Javierre, C. Vuik, F.J. Vermolen, and S. van der Zwaag.  A comparison of numerical models for one-dimensional Stefan problems.  Journal of Computational and Applied Mathematics, 192(2):445–459, 2006.
  * [13] R. Johansson.  Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib.  Apress, 2018.
  * [14] A. Jüngel and A. Unterreiter.  Discrete minimum and maximum principles for finite element approximations of non-monotone elliptic equations.  Numerische Mathematik, 99(3):485–508, 2005.
  * [15] A. Kufner, O. John, and S. Fučík.  Function Spaces.  Springer Science, 1977.
  * [16] K. Kumazaki, T. Aiki, and A. Muntean.  A free boundary problem describing migration into rubbers–quest of the large time behavior.  arXiv preprint arXiv:2102.12766, 2021.
  * [17] K. Kumazaki and A. Muntean.  Local weak solvability of a moving boundary problem describing swelling along a halfline.  Networks & Heterogeneous Media, 14:445–469, 2019.
  * [18] K. Kumazaki and A. Muntean.  Global weak solvability, continuous dependence on data, and large time growth of swelling moving interfaces.  Interfaces and Free Boundaries, 22(1):27–50, 2020.
  * [19] H. G. Landau.  Heat conduction in a melting solid.  Quarterly of Applied Mathematics, 8(1):81–94, 1950.
  * [20] S. Larsson and V. Thomée.  Partial Differential Equations with Numerical Methods, volume 45.  Springer Science & Business Media, 2008.
  * [21] S. Linge and H. P. Langtangen.  Programming for Computations-Python: A Gentle Introduction to Numerical Simulations with Python 3.6.  Springer Nature, 2020.
  * [22] J. A. Mackenzie and M. L. Robertson.  The numerical solution of one-dimensional phase change problems using an adaptive moving mesh method.  Journal of Computational Physics, 161(2):537–557, 2000.
  * [23] M. Morton.  Rubber Technology.  Springer Science & Business Media, 2013.
  * [24] A. Muntean.  Error bounds on semi-discrete finite element approximations of a moving-boundary system arising in concrete corrosion.  International Journal of Numerical Analysis & Modeling, 5(3):353–372, 2008.
  * [25] A. Muntean and M. Böhm.  A moving-boundary problem for concrete carbonation: global existence and uniqueness of weak solutions.  Journal of Mathematical Analysis and Applications, 350(1):234–251, 2009.
  * [26] S. Nepal, R. Meyer, N. H. Kröger, T. Aiki, A. Muntean, Y. Wondmagegne, and U. Giese.  A moving boundary approach of capturing diffusants penetration in to rubber: FEM approximation and comparison with laboratory measurements.  Kautschuk Gummi Kunststoffe, 5:61–69, 2021.
  * [27] J. A. Nitsche.  Finite element approximations to the one dimensional Stefan problem.  In Recent Advances in Numerical Analysis, pages 119–142. Elsevier, 1978.
  * [28] A. K. Pani and P. C. Das.  A finite element method for a single phase semilinear Stefan problem in one space dimension.  Numerical Functional Analysis and Optimization, 12(1-2):153–171, 1991.
  * [29] M. A. Piqueras, R. Company, and L. Jódar.  Numerical analysis and computing of free boundary problems for concrete carbonation chemical corrosion.  Journal of Computational and Applied Mathematics, 336:297–316, 2018\.
  * [30] M. A. Piqueras, R. Company, and L. Jódar.  Solving two-phase freezing Stefan problems: Stability and monotonicity.  Mathematical Methods in the Applied Sciences, 43(14):7948–7960, 2020.
  * [31] M. G. Rezk and J. Foroozesh.  Determination of mass transfer parameters and swelling factor of co2-oil systems at high pressures.  International Journal of Heat and Mass Transfer, 126:380–390, 2018\.
  * [32] V. Thomée.  Galerkin Finite Element Methods for Parabolic Problems, volume 25.  Springer Science & Business Media, 2007.
  * [33] W. Walter.  Differential inequalities and maximum principles: theory, new methods and applications.  Nonlinear Analysis TMA, 8:4695–4711, 1997.
  * [34] E. Zeidler.  Nonlinear Functional Analysis and Its Applications, Linear Monotone Operators, volume 2a.  Springer Verlag, Berlin, 1990.
  * [35] A. Zurek.  Numerical approximation of a concrete carbonation model: Study of the $\sqrt{t}$-law of propagation.  Numerical Methods for Partial Differential Equations, 35(5):1801–1820, 2019.