# GLOBAL $q$-SUPERLINEAR CONVERGENCE OF THE INFINITE-DIMENSIONAL NEWTON’S
METHOD FOR THE REGULARIZED $p$-STOKES EQUATIONS

Niko Schmidt111Christian-Albrecht University of Kiel<EMAIL_ADDRESS>kiel.de

## Abstract

The motion of glaciers can be simulated with the $p$-Stokes equations. We
present an algorithm that solves these equations faster than the Picard
iteration. We do that by proving $q$-superlinear global convergence of the
infinite-dimensional Newton’s method with Armijo step sizes to the solution of
these equations. We only have to add an arbitrarily small diffusion term for
this convergence result. We also consider approximations of exact step sizes.
Exact step sizes are possible because we reformulate the problem as minimizing
a convex functional. Next, we prove that the additional diffusion term only
causes minor differences in the solution compared to the original $p$-Stokes
equations. Finally, we test our algorithms on a reformulation of the
experiment ISMIP-HOM $B$. The approximation of exact step sizes for the Picard
iteration and Newton’s method is superior in the experiment compared to the
Picard iteration. Also, Newton’s method with Armijo step sizes converges
faster than the Picard iteration. However, the reached accuracy of Newton’s
method with Armijo step sizes depends more on the resolution of the domain.

### Keywords

$p$-Stokes, Newton’s method, global convergence, glaciology

##  1\. Introduction

Ice-sheet models spend most of the computation time solving the momentum
equations, [1]. These equations are nonlinear partial differential equations
named $p$-Stokes equations.

We prove that Newton’s method with Armijo step sizes is globally
$q$-superlinear convergent to the solution of the $p$-Stokes equations if we
add a small diffusion term. We control the step size with a convex functional,
which is the anti-derivative of the $p$-Stokes equations. Evaluating this
functional has nearly no computational cost. Moreover, the convexity of the
functional allows us to approximate the exact step sizes. Thus, it seems
superior to approaches that use the residual norm for step size control. We
slightly modify the numerical experiment ISMIP-HOM $B$ to test Newton’s method
with Armijo step sizes. Additionally, we test the Picard iteration and
Newton’s method with approximations of exact step sizes. The small diffusion
term is necessary to verify Fréchet differentiability of the $p$-Stokes
equations in the variational formulation. Furthermore, the shear-thinning
viscosity term has a negative exponent. Thus, we need for differentiability a
positive constant in this term. To conclude the theory, we show that the
regularized solution converges to the solution of the $p$-Stokes equations in
a Lebesgue space for vanishing regularization.

In the beginning, we need to verify the well-posedness of the equations. There
are different types of literature results for well-posedness: complex boundary
conditions [2], an implicitly given viscosity [3] with a differentiable shear-
thinning term, or a differentiable shear-thinning term with Dirichlet boundary
conditions, [4]. One recent publication with more general boundary conditions
than we consider in two dimensions uses Newton’s method in a finite-
dimensional setting, see [5]. However, we need a combination of a
differentiable, explicitly given shear-thinning viscosity and complex boundary
conditions. Glaciologists use such a formulation to simulate glaciers, [6].
Gâteaux differentiability results are mainly devoted to optimal control and
the necessary control-to-state mapping, see for example, [7] for the
$p$-Laplace equations or [8] for the Navier-Stokes equations with a nonlinear
$p$-Stokes term. In [8], the Gâteaux differentiability is shown for $p\geq 2$.
The idea of how differentiability can be shown for $1<p<2$ is briefly
mentioned in [9, section $6$]. This paper uses the additional diffusion term
for optimal control. It proves convergence in Sobolev spaces for vanishing
regularization of the diffusion term. However, that paper has different
requirements: It considers optimal control in a slightly different
formulation, has Dirichlet boundary conditions, and has a convective term. It
has more restrictions on the exponent of the shear-thinning fluids than we
have. In [3], local quadratic convergence of Newton’s method is shown for the
finite-dimensional case. Glaciologists already consider Newton’s method, [10],
but we consider a different approach by adding a small diffusive term and
using a convex functional for step size control. There is also a multigrid
approach with Newton’s method and a reformulation of the partial differential
equation with first derivatives available, [11].

The paper is structured as follows: In section 2, we derive the variational
formulation, project to divergence-free spaces, and verify existence and
uniqueness of a solution. In section 3, we consider Fréchet differentiability
and introduce the regularized $p$-Stokes equations. In section 4 and section
5, we verify local and global $q$-superlinear convergence, respectively. For
the latter, we formulate the problem as a minimization problem. In section 6,
we verify that the solution of the regularized $p$-Stokes equations converges
to the solution of the $p$-Stokes equations for vanishing regularization. In
section 7, we consider a numerical experiment. We summarize our results in the
final section 8.

##  2\. The p-Stokes equations

In this section, we formulate the $p$-Stokes equations in both the classical
and the variational formulations. The $p$-Stokes equations model the motion of
glaciers. Their complexity results from nonlinear viscosity, also described as
shear-thinning, [6]. Mass conservation and incompressibility lead to
divergence-free velocity $\boldsymbol{v}$. Let $N\in\\{2,3\\}$,
$\Omega\subseteq\mathbb{R}^{N}$ be a Lipschitz domain. Let us formally write
down the $p$-Stokes equations without boundary conditions now. Let $\sigma$ be
the stress tensor, $\boldsymbol{v}$ the velocity, $\rho$ the density, and
$\boldsymbol{g}$ the gravitational acceleration. The $p$-Stokes equations are:

$\displaystyle-\mathrm{div}\sigma$
$\displaystyle=-\rho\boldsymbol{g}\quad\text{on }\Omega,$
$\displaystyle\mathrm{div}\boldsymbol{v}$ $\displaystyle=0\quad\quad\text{ on
}\Omega$

by using $\sigma$ as a second-order tensor and the divergence on it column-
wise. The relation between the velocity $\boldsymbol{v}$, the pressure $\pi$,
and the stress tensor $\sigma$ is given by

$\displaystyle\sigma=-\pi I+BS^{p}(D\boldsymbol{v})\text{,}\quad B\in
L^{\infty}(\Omega)\text{, }\quad B\geq c_{0}\in(0,\infty)$

with the identity matrix $I=(I_{ij})_{ij}$, the mapping
$S^{p}:\mathbb{R}^{N\times N}\to\mathbb{R}^{N\times N}$,

$\displaystyle S^{p}(P)$
$\displaystyle=\big{(}|P|^{2}+\delta^{2}\big{)}^{(p-2)/2}P\text{, }\quad
P\in\mathbb{R}^{N\times N}\text{, }\quad\delta>0\text{, }\quad p\in(1,2)$

with the Frobenius norm $|\cdot|$ and

$\displaystyle(D\boldsymbol{v})_{ij}$
$\displaystyle=\frac{1}{2}\bigg{(}\frac{\partial v_{i}}{\partial
x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\bigg{)}\text{, }\quad
i,j\in\\{1,...,N\\}.$

We interpret $S^{p}$ for vectors $\boldsymbol{v}\in\mathbb{R}^{N}$ as a
mapping from $\mathbb{R}^{N}$ to $\mathbb{R}^{N}$.

The case $p=2$ reduces to the Stokes problem. In glaciological applications
$p=4/3$, see [6, section 1.4.3], or in more recent approaches $p=5/4$, see
[12, Abstract], is used. The following result states information about the
integrability of $S^{p}$:

######  Lemma 2.1.

For all $Q\in L^{p}(\Omega)^{N\times N}$ follows $S^{p}(Q)\in
L^{q}(\Omega)^{N\times N}$ with the dual exponent $q$.

###### Proof.

With the dual exponent $q=p/(p-1)$ follows

$\displaystyle|S^{p}(Q)|\leq|Q|^{p-1}\Rightarrow|S^{p}(Q)|^{q}\leq|Q|^{p}\in
L^{1}(\Omega).$

∎

The Dirichlet boundary condition $\boldsymbol{v}=0$ is satisfied for those
parts of the glacier frozen to the ground. The interaction of the glacier with
the air is given by $\sigma\cdot\boldsymbol{n}=0$ with the matrix-vector
multiplication of the matrix $\sigma$ and the unit normal vector
$\boldsymbol{n}$. We assume nonlinear sliding at parts of the bedrock that are
not frozen to the ground. This sliding is represented by tangential sliding.
For these parts of the bedrock, we assume that the normal component of the
velocity is $0$ because we neglect the melting of ice or freezing of water at
the bedrock. The boundary conditions are as follows:

$\displaystyle\boldsymbol{v}$
$\displaystyle=0\quad\quad\quad\quad\quad\text{on }\Gamma_{d},$
$\displaystyle\sigma\cdot\boldsymbol{n}$
$\displaystyle=0\quad\quad\quad\quad\quad\text{on }\Gamma_{a}$
$\displaystyle\boldsymbol{\sigma}_{t}$ $\displaystyle=-\tau
S^{s}(\boldsymbol{v_{t}})\quad\text{ on }\Gamma_{b}\text{,}\quad\tau\in
L^{\infty}(\Gamma_{b})\text{, }\quad\tau\geq 0\text{, }\quad s\in(1,p],$
$\displaystyle\boldsymbol{v}\cdot\boldsymbol{n}$
$\displaystyle=0\quad\quad\quad\quad\quad\text{on }\Gamma_{b},$

with $\partial\Omega=\Gamma_{d}\cup\Gamma_{a}\cup\Gamma_{b}$,
$|\Gamma_{d}|>0$, the definition of $\sigma:=-\pi I+BS^{p}(D\boldsymbol{v})$,
and the relation between tangential and normal components:

(1)
$\displaystyle\boldsymbol{v}=\boldsymbol{v_{t}}+(\boldsymbol{v}\cdot\boldsymbol{n})\boldsymbol{n}\quad\text{
and
}\quad\sigma\cdot\boldsymbol{n}=\boldsymbol{\sigma_{t}}+(\sigma\cdot\boldsymbol{n}\cdot\boldsymbol{n})\boldsymbol{n}.$

The restriction $|\Gamma_{d}|>0$ is also fulfilled in common applications, see
for example [12, Abstract] or [13, Abstract]. The boundary conditions imply as
a natural choice of the space

$\displaystyle W_{r}:=\\{\boldsymbol{v}\in W^{1,r}(\Omega)^{N}\text{;
}\boldsymbol{v}|_{\Gamma_{d}}=0\text{,
}\boldsymbol{v}|_{\Gamma_{b}}\cdot\boldsymbol{n}|_{\Gamma_{b}}=0\\}\text{,
}\quad r>1$

with norm

$\displaystyle\|\boldsymbol{v}\|_{W_{r}}^{r}=\int_{\Omega}\bigg{(}\sum_{i,j=1}^{N}\Big{(}\frac{\partial
v_{i}}{\partial x_{j}}\Big{)}^{2}\bigg{)}^{r/2}\,dx.$

The Poincaré inequality implies that the $W_{r}$-norm is equivalent to the
$W^{1,r}(\Omega)^{N}$-norm. For example, the Poincaré inequality was proved in
[14, Theorem 1.5] for $p=2$. But the proof is identical for $p\in(1,2)$.

###  2.1. Variational formulation

We derive the variational formulation in the space $W_{r}$. We define the
double scalar product and obtain the following equation for
$\sigma,\tau\in\mathbb{R}^{N\times N}$ immediately:

(2) $\displaystyle\sigma:\tau:=\sum_{i,j=1}^{N}\sigma_{ij}\tau_{ij}\text{,
}\quad\sigma:\tau^{T}=\sum_{i,j=1}^{N}\sigma_{ij}\tau_{ji}=\sigma^{T}:\tau.$

Let $r^{\prime}$ be the dual exponent, $\boldsymbol{\phi}\in W_{r^{\prime}}$,
$\sigma\in W^{1,r}(\Omega)^{N\times N}$. We conclude with partial integration

(3)
$\displaystyle-\int_{\Omega}\mathrm{div}\sigma\cdot\boldsymbol{\phi}\,dx=\int_{\Omega}\sigma:\nabla\boldsymbol{\phi}\,dx+\int_{\partial\Omega}\sigma\cdot\boldsymbol{n}\cdot\boldsymbol{\phi}\,ds.$

Let $\boldsymbol{v}\in W_{r}$, $\pi\in L^{r^{\prime}}(\Omega)$. We use
$\sigma=-\pi I+BS^{p}(D\boldsymbol{v})$ and
$I:\nabla\boldsymbol{\phi}=\mathrm{div}\boldsymbol{\phi}$ to obtain

$\displaystyle\int_{\Omega}\sigma:\nabla\boldsymbol{\phi}\,dx=-\int_{\Omega}\pi\mathrm{div}\boldsymbol{\phi}\,dx+\int_{\Omega}BS^{p}(D\boldsymbol{v}):\nabla\boldsymbol{\phi}\,dx$

for the first summand on the right-hand side of relation
$(\ref{VarSummands})$. The second summand on the right-hand side vanishes on
$\Gamma_{d}$ because the test function $\boldsymbol{\phi}\in W_{r^{\prime}}$
vanishes on $\Gamma_{d}$. Moreover, $\sigma\cdot\boldsymbol{n}=0$ is valid on
$\Gamma_{a}$. Thus, the integral over this domain disappears, too. The
relation between tangential and normal components, see $(\ref{Tangential})$
implies

$\displaystyle(\sigma\cdot\boldsymbol{n})\cdot\boldsymbol{\phi}=(\sigma\cdot\boldsymbol{n}\cdot\boldsymbol{n})(\boldsymbol{\phi}\cdot\boldsymbol{n})+\boldsymbol{\sigma_{t}}\cdot\boldsymbol{\phi_{t}}.$

We split the boundary integral on $\Gamma_{b}$ into normal and tangential
components and use the boundary conditions on $\Gamma_{b}$ to obtain

$\displaystyle-\int_{\Gamma_{b}}\sigma\cdot\boldsymbol{n}\cdot\boldsymbol{\phi}\,ds$
$\displaystyle=-\int_{\Gamma_{b}}(\sigma\cdot\boldsymbol{n}\cdot\boldsymbol{n})(\boldsymbol{\phi}\cdot\boldsymbol{n})\,ds-\int_{\Gamma_{b}}\boldsymbol{\sigma}_{t}\cdot\boldsymbol{\phi}_{t}\,ds$
$\displaystyle=\int_{\Gamma_{b}}\tau\big{(}|\boldsymbol{v_{t}}|^{2}+\delta^{2}\big{)}^{(s-2)/2}\boldsymbol{v_{t}}\cdot\boldsymbol{\phi}_{t}\,ds.$

Due to equation $(\ref{Tangential})$, we have
$\boldsymbol{v}=\boldsymbol{v_{t}}$ on $\Gamma_{b}$, which yields

$\displaystyle\int_{\Gamma_{b}}\tau\big{(}|\boldsymbol{v_{t}}|^{2}+\delta^{2}\big{)}^{(s-2)/2}\boldsymbol{v_{t}}\cdot\boldsymbol{\phi}_{t}\,ds$
$\displaystyle=\int_{\Gamma_{b}}\tau\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(s-2)/2}\boldsymbol{v}\cdot\boldsymbol{\phi}\,ds$
$\displaystyle=\int_{\Gamma_{b}}\tau
S^{s}(\boldsymbol{v})\cdot\boldsymbol{\phi}\,ds.$

We set $L^{q}_{0}(\Omega):=\\{\pi\in L^{q}(\Omega)\text{;
}\int_{\Omega}\pi\,dx=0\\}$. Then a weak solution of the problem is given by
$(\boldsymbol{v},\pi)\in(W_{p},L^{q}_{0}(\Omega))$ with

$\displaystyle\langle
A\boldsymbol{v},\boldsymbol{\phi}\rangle_{W_{p}^{*},W_{p}}-(\pi,\mathrm{div}\boldsymbol{\phi})$
$\displaystyle=-(\rho\boldsymbol{g},\boldsymbol{\phi})\quad\text{for all
}\boldsymbol{\phi}\in W_{p},$
$\displaystyle-(\mathrm{div}\boldsymbol{v},\psi)$
$\displaystyle=0\qquad\qquad\text{ for all }\psi\in L^{q}_{0}(\Omega)$

with the operator $A:W_{p}\to W_{p}^{*}$,

$\displaystyle\langle
A\boldsymbol{v},\boldsymbol{\phi}\rangle_{V_{p}^{*},V_{p}}=\int_{\Omega}BS^{p}(D\boldsymbol{v}):\nabla\boldsymbol{\phi}\,dx+\int_{\Gamma_{b}}\tau
S^{s}(\boldsymbol{v})\cdot\boldsymbol{\phi}\,ds\text{,
}\quad\boldsymbol{\phi}\in W_{p}.$

The operator $A$ is well-defined due to Lemma 2.1. We set
$V_{p}:=\\{\boldsymbol{v}\in W_{p}\text{; }\mathrm{div}\boldsymbol{v}=0\\}$ to
consider divergence-free velocities. Now, we examine weak solutions of
$A:V_{p}\to V_{p}^{*}$,

(4) $\displaystyle\langle
A\boldsymbol{v},\boldsymbol{\phi}\rangle_{V_{p}^{*},V_{p}}=-\int_{\Omega}\rho\boldsymbol{g}\cdot\boldsymbol{\phi}\,dx\quad\text{
for all }\phi\in V_{p}.$

By solving the divergence-free problem, we solve the original problem: The
inf-sup condition is proved by [15, Corollary 3.2] for Dirichlet boundary
conditions:

$\displaystyle\inf_{\pi\in L^{q}_{0}(\Omega)}\sup_{\boldsymbol{v}\in
W^{1,p}_{0}(\Omega)^{N}}\frac{(\pi,\mathrm{div}\boldsymbol{v})}{\|\boldsymbol{v}\|_{W^{1,p}_{0}(\Omega)^{N}}\|\pi\|_{L^{q}(\Omega)}}\geq
c>0.$

The inf-sup condition is also fulfilled for $W_{p}\supseteq
W^{1,p}_{0}(\Omega)^{N}$. Thus, we can apply [16, Theorem IV.1.4]. This yields
for each solution in the divergence-free formulation a unique solution for the
original problem. This Theorem is stated in Hilbert spaces, but the proof is
identical for operators on Banach spaces with an existing dual operator.

###  2.2. Existence and uniqueness of weak solutions in the divergence-free
space

In this subsection, we verify conditions to conclude existence and uniqueness
of solutions. These conditions are formulated in the following theorem.

######  Theorem 2.2 (Browder-Minty Theorem).

Let $X$ be a reflexive separable Banach space and $A:X\to X^{*}$ a strictly
monotone, coercive, and continuous operator. Let $f\in X^{*}$. Then there
exists a unique solution $u\in X$ for the equation

$\displaystyle Au=f\quad\text{ in }X^{*}.$

###### Proof.

The existence of $u\in X$ is proved in [17] Theorem $2$ in a more general
version. Uniqueness follows immediately with the strict monotonicity of $A$. ∎

We analyze $S^{p}$ and $S^{s}$ to verify that $A$ is strictly monotone,
coercive, and continuous. For that purpose, we need the following result:

######  Lemma 2.3.

Let $S:\mathbb{R}^{N\times N}\to\mathbb{R}^{N\times N}$ with

(5) $\displaystyle\sum_{i,j,k,\ell=1}^{N}\frac{\partial S_{ij}(P)}{\partial
P_{kl}}Q_{ij}Q_{k\ell}$ $\displaystyle\geq c_{1}(|P|+\delta)^{p-2}|Q|^{2},$
(6) $\displaystyle\bigg{|}\frac{\partial S_{ij}(P)}{\partial
P_{k\ell}}\bigg{|}$ $\displaystyle\leq c_{2}(|P|+\delta)^{p-2}\quad\text{ for
all }i,j,k,\ell\in\\{1,...,N\\}$

for all $P,Q\in\mathbb{R}^{N\times N}$, $c_{1},c_{2}\in(0,\infty)$,
$\delta\geq 0$ and $p>1$. Then, there exist $c,C\in\mathbb{R}$ independent of
$p$ for all $P,Q\in\mathbb{R}^{N\times N}$ with

$\displaystyle(S(P)-S(Q)):(P-Q)$ $\displaystyle\geq
c(\delta+|P|+|Q|)^{p-2}|P-Q|^{2},$ $\displaystyle|S(P)-S(Q)|$
$\displaystyle\leq C(\delta+|P|+|Q|)^{p-2}|P-Q|.$

###### Proof.

See [18, Lemma 6.3]. ∎

We calculate the derivative of $S^{r}$, $r\in\\{s,p\\}$ and verify the
properties stated above to apply the Browder-Minty Theorem 2.2:

######  Lemma 2.4.

Let $r\in\\{s,p\\}$. The function $S^{r}_{kl}$ is continuously differentiable.
Furthermore, the inequalities $(\ref{LowerBound})$ and $(\ref{UpperBound})$
are fulfilled for $S^{r}$.

###### Proof.

We verify the conditions in Lemma 2.3. Let $P\in\mathbb{R}^{N\times N}$,
$i,j\in\\{1,...,N\\}$. We calculate

$\displaystyle\frac{\partial}{\partial
P_{ij}}\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}=(r-2)\Bigg{(}\bigg{(}\sum_{k,\ell=1}^{N}P_{k\ell}^{2}\bigg{)}+\delta^{2}\Bigg{)}^{(r-4)/2}P_{ij}.$

Let $k,\ell\in\\{1,..,N\\}$. We infer

(7) $\displaystyle\begin{split}\frac{\partial S^{r}_{k\ell}}{\partial
P_{ij}}(P)&=\frac{\partial}{\partial
P_{ij}}\Big{(}\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}P_{k\ell}\Big{)}\\\
&=(r-2)\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-4)/2}P_{ij}P_{k\ell}+\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}I_{ik}I_{j\ell}.\end{split}$

We can verify inequality $(\ref{UpperBound})$ with
$\tilde{c}_{1}\in\mathbb{R}$ now:

$\displaystyle\bigg{|}\frac{\partial S^{r}_{k\ell}(P)}{\partial
P_{ij}}\bigg{|}$ $\displaystyle\leq
N^{2}\Big{(}(2-r)\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-4)/2}|P|^{2}+\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\Big{)}$
$\displaystyle=N^{2}\Big{(}(2-r)\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\frac{|P|\,|P|}{|P|^{2}+\delta^{2}}+\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\Big{)}$
$\displaystyle\leq\tilde{c}_{1}N^{2}(3-r)(|P|+\delta)^{r-2}.$

We verify inequality $(\ref{LowerBound})$ next. Let $Q\in\mathbb{R}^{N\times
N}$. We conclude

$\displaystyle\quad\sum_{i,j,k,\ell=1}^{N}\frac{\partial}{\partial
P_{ij}}\big{(}S^{r}_{k\ell}(P)\big{)}Q_{k\ell}Q_{ij}$
$\displaystyle=\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\bigg{(}\sum_{i,j,k,\ell=1}^{N}(r-2)\frac{P_{ij}P_{k\ell}}{|P|^{2}+\delta^{2}}Q_{k\ell}Q_{ij}+\sum_{i,j,k,\ell=1}^{N}I_{ik}I_{jl}Q_{k\ell}Q_{ij}\bigg{)}$
$\displaystyle=\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\bigg{(}(r-2)\frac{1}{|P|^{2}+\delta^{2}}\sum_{k,\ell=1}^{N}P_{k\ell}Q_{k\ell}\sum_{i,j=1}^{N}P_{ij}Q_{ij}+\sum_{i,j=1}^{N}Q_{ij}^{2}\bigg{)}$
$\displaystyle=\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\bigg{(}(r-2)\frac{(P:Q)^{2}}{|P|^{2}+\delta^{2}}+|Q|^{2}\bigg{)}.$

We conclude with $r-2\leq 0$, $(P:Q)^{2}\leq|P|^{2}|Q|^{2}$, and
$\tilde{c}_{2}\in\mathbb{R}$

$\displaystyle\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}\bigg{(}(r-2)\frac{(P:Q)^{2}}{|P|^{2}+\delta^{2}}+|Q|^{2}\bigg{)}$
$\displaystyle\geq\big{(}|P|^{2}+\delta^{2}\big{)}^{(r-2)/2}(r-1)|Q|^{2}$
$\displaystyle\geq\tilde{c}_{2}(r-1)(|P|+\delta)^{r-2}|Q|^{2}.$

We set $P:=I_{ij}p$ and $Q:=I_{ij}q$ with $p,q\in\mathbb{R}^{N}$ for the
vector-valued situation. ∎

We verify more than the continuity of the operator $A:V_{p}\to V_{p}^{*}$:

######  Lemma 2.5.

The operator $A:V_{p}\to V_{p}^{*}$ is Hölder continuous with exponent $p-1$.

###### Proof.

Let $\boldsymbol{v},\boldsymbol{w},\boldsymbol{\phi}\in V_{p}$. We obtain for
the second summand of the operator $A$ with Lemma 2.4 and $C\in\mathbb{R}$

$\displaystyle\bigg{|}\int_{\Gamma_{b}}\tau
S^{s}(\boldsymbol{v}-\boldsymbol{w})\cdot\boldsymbol{\phi}\,ds\bigg{|}$
$\displaystyle\leq
C\int_{\Gamma_{b}}(\delta+|\boldsymbol{v}|+|\boldsymbol{w}|)^{s-2}|\boldsymbol{v}-\boldsymbol{w}|\,|\boldsymbol{\phi}|\,ds$
$\displaystyle=C\int_{\Gamma_{b}}(\delta+|\boldsymbol{v}|+|\boldsymbol{w}|)^{s-p}\frac{|\boldsymbol{v}-\boldsymbol{w}|}{(\delta+|\boldsymbol{v}|+|\boldsymbol{w}|)^{2-p}}|\boldsymbol{\phi}|\,ds$
$\displaystyle\leq
C\int_{\Gamma_{b}}\qquad\delta^{s-p}\qquad\qquad\quad\frac{|\boldsymbol{v}-\boldsymbol{w}|}{|\boldsymbol{v}-\boldsymbol{w}|^{2-p}}\qquad|\boldsymbol{\phi}|\,ds$
$\displaystyle=C\delta^{s-p}\int_{\Gamma_{b}}|\boldsymbol{v}-\boldsymbol{w}|^{p-1}|\boldsymbol{\phi}|\,ds$
$\displaystyle\leq
C\delta^{s-p}\bigg{(}\int_{\Gamma_{b}}|\boldsymbol{v}-\boldsymbol{w}|^{q(p-1)}\,ds\bigg{)}^{1/q}\|\boldsymbol{\phi}\|_{L^{p}(\Gamma_{b})}.$

Thus, we obtain with the dual exponent $q=p/(p-1)$

$\displaystyle\bigg{(}\int_{\Gamma_{b}}|\boldsymbol{v}-\boldsymbol{w}|^{q(p-1)}\,ds\bigg{)}^{1/q}=\bigg{(}\int_{\Gamma_{b}}|\boldsymbol{v}-\boldsymbol{w}|^{p}\,ds\bigg{)}^{(p-1)/p}=\|\boldsymbol{v}-\boldsymbol{w}\|_{L^{p}(\Gamma_{b})}^{p-1}.$

We infer with the trace operator, [19, Theorem 1.12], and $c\in\mathbb{R}$

$\displaystyle\|\boldsymbol{v}-\boldsymbol{w}\|_{L^{p}(\Gamma_{b})}^{p-1}\|\boldsymbol{\phi}\|_{L^{p}(\Gamma_{b})}\leq
c\|\boldsymbol{v}-\boldsymbol{w}\|_{V_{p}}^{p-1}\|\boldsymbol{\phi}\|_{V_{p}}.$

This upper bound yields the Hölder continuity with exponent $p-1$. The same
arguments are valid for the first summand of the operator $A$ with $\Omega$
instead of $\Gamma_{b}$, $B$ instead of $\tau$, and $s=p$. ∎

We need the following result to prove coercivity:

######  Lemma 2.6 (Korn’s inequality).

For $\boldsymbol{v}\in V_{p}$ we have with $c\in\mathbb{R}$

$\displaystyle\|\boldsymbol{v}\|_{V_{p}}\ \leq
c\|D\boldsymbol{v}\|_{L^{p}(\Omega)^{N\times N}}.$

###### Proof.

See [2, Appendix Lemma 3] for $N=3$. The case $N=2$ follows with the constant
continuation of $\boldsymbol{v}=(v_{1},v_{2})$ to the third dimension. ∎

We obtain a lower bound in the following sense:

######  Lemma 2.7.

The operator $A$ is coercive.

###### Proof.

We conclude for the second summand of the operator $A$

$\displaystyle\int_{\Gamma_{b}}\tau\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(s-2)/2}\boldsymbol{v}\cdot\boldsymbol{v}\,ds\geq
0.$

Equation $(\ref{Transpose})$ gives
$D\boldsymbol{v}:\nabla\boldsymbol{v}=D\boldsymbol{v}:D\boldsymbol{v}=|D\boldsymbol{v}|^{2}$.
In [20, Lemma 1.19 in chapter 5], it is shown
$(|D\boldsymbol{v}|^{2}+1)^{(p-2)/2}|D\boldsymbol{v}|^{2}\geq
c|D\boldsymbol{v}|^{p}$ with $c\in\mathbb{R}$. Thus, we have:

$\displaystyle\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}|D\boldsymbol{v}|^{2}\,dx$
$\displaystyle=\int_{\Omega}B\delta^{p}\bigg{(}\frac{|D\boldsymbol{v}|^{2}}{\delta^{2}}+1\bigg{)}^{(p-2)/2}\bigg{|}\frac{D\boldsymbol{v}}{\delta}\bigg{|}^{2}\,dx$
$\displaystyle\geq\int_{\Omega}Bc\delta^{p}\bigg{|}\frac{D\boldsymbol{v}}{\delta}\bigg{|}^{p}\,dx=\int_{\Omega}Bc|D\boldsymbol{v}|^{p}\,dx.$

Korn’s inequality and $B\geq c_{0}$ yield $\tilde{c}\in\mathbb{R}$ with

(8)
$\displaystyle\frac{\displaystyle\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}|D\boldsymbol{v}|^{2}\,dx}{\|\boldsymbol{v}\|_{V_{p}}}\geq\frac{\tilde{c}\|\boldsymbol{v}\|_{V_{p}}^{p}}{\|\boldsymbol{v}\|_{V_{p}}}=\tilde{c}\|\boldsymbol{v}\|_{V_{p}}^{p-1}\to\infty\text{
for }\|\boldsymbol{v}\|_{V_{p}}\to\infty.$

∎

Hence, we can prove our existence result:

######  Theorem 2.8.

There exists exactly one weak solution of the divergence-free $p$-Stokes
equations, see equation $(\ref{WeakSolution})$.

###### Proof.

(In [4] the case $\Gamma_{d}=\partial\Omega$ is analyzed.) Lemma 2.5 yields
that the operator $A$ is Hölder continuous and thus continuous. We verified
the coercivity of $A$ in Lemma 2.7. The strict monotonicity of $S^{p}$ and
$S^{s}$, proved in Lemma 2.4, yield the strict monotonicity of $A$: Let
$\boldsymbol{v},\boldsymbol{w}\in V_{p}$. We have

$\displaystyle\langle
A\boldsymbol{v}-A\boldsymbol{w},\boldsymbol{v}-\boldsymbol{w}\rangle_{V_{p}^{*},V_{p}}$
$\displaystyle=$
$\displaystyle\int_{\Omega}B\big{(}S^{p}(D\boldsymbol{v})-S^{p}(D\boldsymbol{w})\big{)}:(\nabla\boldsymbol{v}-\nabla\boldsymbol{w})\,dx+\int_{\Gamma_{b}}\tau\big{(}S^{s}(\boldsymbol{v})-S^{s}(\boldsymbol{w})\big{)}\cdot(\boldsymbol{v}-\boldsymbol{w})\,ds$
$\displaystyle\geq$
$\displaystyle\int_{\Omega}Bc(\delta+|D\boldsymbol{v}|+|D\boldsymbol{w}|)^{p-2}|D\boldsymbol{v}-D\boldsymbol{w}|^{2}\,dx+\int_{\Gamma_{b}}\tau(\delta+|\boldsymbol{v}|+|\boldsymbol{w}|)^{p-2}|\boldsymbol{v}-\boldsymbol{w}|\,ds$
$\displaystyle\geq$ $\displaystyle
c_{0}c\int_{\Omega}|D\boldsymbol{v}-D\boldsymbol{w}|^{p}\,dx$
$\displaystyle\geq$ $\displaystyle
c_{0}\tilde{c}\|\boldsymbol{v}-\boldsymbol{w}\|_{V_{p}}^{p}>0.$

The last inequality followed with Korn’s inequality, see Lemma 2.6 and
$\tilde{c}\in\mathbb{R}$. We conclude the existence and uniqueness of a
solution with the Browder-Minty Theorem, see 2.2. ∎

##  3\. Differentiability

We must calculate the Fréchet derivative of the operator $A$ to apply Newton’s
method. At first, we only consider the first summand of the operator $A$.

###  3.1. Fréchet differentiability on a subset

We show differentiability in the space $V_{2}$. We will see that this
derivative is not well-posed in $V_{p}$. Let $G_{1}:V_{2}\to V_{2}^{*}$,

$\displaystyle\langle
G_{1}\boldsymbol{v},\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}:=(BS^{p}(D\boldsymbol{v}),\nabla\boldsymbol{\phi})\text{,
}\quad\boldsymbol{\phi}\in V_{2}.$

The operator $G_{1}$ has a directional derivative:

######  Lemma 3.1 (Directional derivative).

For all $\boldsymbol{v}\in V_{2}$ the function $G_{1}$ is differentiable in
direction $\boldsymbol{w}\in V_{2}$ with

(9) $\displaystyle\begin{split}\langle
G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}&=\int_{\Omega}(p-2)B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-4)/2}(D\boldsymbol{v}:D\boldsymbol{w})\,(D\boldsymbol{v}:\nabla\boldsymbol{\phi})\,dx\\\
&\quad+\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}D\boldsymbol{w}:\nabla\boldsymbol{\phi}\,dx\text{,
}\quad\boldsymbol{\phi}\in V_{2}.\end{split}$

###### Proof.

The $o$-notation is the limit to zero in the following calculations. We prove
pointwise convergence. For all $i,j\in\\{1,...,N\\}$, we calculate the Taylor
expansion of $S^{p}$ in $D\boldsymbol{v}$ with the continuous derivative, see
equation $(\ref{DerivativeSr})$. This yields almost everywhere

$\displaystyle S^{p}_{ij}(D\boldsymbol{v}+tD\boldsymbol{w})$
$\displaystyle=S^{p}_{ij}(D\boldsymbol{v})+\sum_{k,\ell=1}^{N}(p-2)\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-4)/2}(D\boldsymbol{v})_{ij}(D\boldsymbol{v})_{k\ell}t(D\boldsymbol{w})_{k\ell}$
$\displaystyle\quad+\sum_{k,\ell=1}^{N}\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}I_{ik}I_{j\ell}t(D\boldsymbol{w})_{k\ell}+o(|tD\boldsymbol{w}|)$
$\displaystyle=S^{p}_{ij}(D\boldsymbol{v})+t(p-2)\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-4)/2}(D\boldsymbol{v})_{ij}D\boldsymbol{v}:D\boldsymbol{w}$
$\displaystyle\quad+(|D\boldsymbol{v}|^{2}+\delta^{2})^{(p-2)/2}t(D\boldsymbol{w})_{ij}+o(|tD\boldsymbol{w}|).$

This implies for $t\to 0$

$\displaystyle\quad\frac{B}{t}\big{(}S^{p}_{ij}(D\boldsymbol{v}+tD\boldsymbol{w})-S^{p}_{ij}(D\boldsymbol{v})\big{)}$
$\displaystyle=(p-2)B(|D\boldsymbol{v}|^{2}+\delta^{2})^{(p-4)/2}(D\boldsymbol{v})_{ij}D\boldsymbol{v}:D\boldsymbol{w}+B(|D\boldsymbol{v}|^{2}+\delta^{2})^{(p-2)/2}(D\boldsymbol{w})_{ij}$

for all $i,j\in\\{1,...,N\\}$ almost everywhere. We obtained the pointwise
convergence. To apply the dominated convergence theorem, we calculate a
majorant. Using the Lipschitz continuity of $S^{p}$, see Lemma 2.4, and
$C\in\mathbb{R}$, we conclude

$\displaystyle\quad\bigg{|}\int_{\Omega}\frac{B}{t}\Big{(}S^{p}(D\big{(}\boldsymbol{v}+t\boldsymbol{w}\big{)}\big{)}-S^{p}(D\boldsymbol{v})\Big{)}:\nabla\phi\,dx\bigg{|}$
$\displaystyle\leq
C\int_{\Omega}\frac{B}{t}(\delta+|D(\boldsymbol{v}+t\boldsymbol{w})|+|D\boldsymbol{v}|)^{p-2}|D(\boldsymbol{v}+t\boldsymbol{w})-D\boldsymbol{v}|\,|\nabla\boldsymbol{\phi}|\,dx$
$\displaystyle\leq
C\delta^{p-2}\|B\|_{L^{\infty}(\Omega)}\int_{\Omega}|D\boldsymbol{w}|\,|\nabla\boldsymbol{\phi}|\,dx$
(10) $\displaystyle\leq
C\delta^{p-2}\|B\|_{L^{\infty}(\Omega)}\|\boldsymbol{w}\|_{V_{2}}\|\boldsymbol{\phi}\|_{V_{2}}<\infty.$

∎

The operator $G_{1}$ has only a directional derivative on $V_{2}$ not on
$V_{p}$:

######  Remark 3.2.

The second summand of $G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w})$ in Lemma
3.1 is not well-defined for all $\boldsymbol{w},\boldsymbol{\phi}\in V_{p}$.

###### Proof.

We know $D\boldsymbol{w},\nabla\boldsymbol{\phi}\in L^{p}(\Omega)^{N\times N}$
for all $\boldsymbol{w},\boldsymbol{\phi}\in V_{p}$. Thus, the integral is not
defined for example for $B\equiv 1$, $\boldsymbol{v}\equiv 0$, and suitable
$\boldsymbol{w},\boldsymbol{\phi}\in V_{p}$. ∎

We verify Fréchet differentiability of the restricted problem now:

######  Theorem 3.3.

The operator $G_{1}:V_{2}\to V_{2}^{*}$ is Fréchet differentiable with the
derivative as in equation $(\ref{DerivativeFormula})$. We set
$G_{1}^{\prime}(\boldsymbol{v})\boldsymbol{w}:=G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w})$
to emphasize the Fréchet differentiability.

###### Proof.

We calculated the directional derivative in Lemma 3.1. Obviously, the operator
$\boldsymbol{w}\mapsto\langle
G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}$
is linear. We easily obtain the boundedness with relation $(\ref{Bound})$ in
Lemma 3.1:

$\displaystyle|\langle
G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}|$
$\displaystyle=\bigg{|}\lim_{t\downarrow
0}\frac{G_{1}(\boldsymbol{v}+t\boldsymbol{w})-G_{1}(\boldsymbol{v}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}}{t}\bigg{|}$
$\displaystyle=\bigg{|}\lim_{t\downarrow
0}\int_{\Omega}\frac{B}{t}\big{(}S^{p}(D(\boldsymbol{v}+t\boldsymbol{w}))-S^{p}(D\boldsymbol{v})\big{)}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
(11) $\displaystyle\leq
C\delta^{p-2}\|B\|_{L^{\infty}(\Omega)}\|\boldsymbol{w}\|_{V_{2}}\|\boldsymbol{\phi}\|_{V_{2}}.$

Thus, $\boldsymbol{w}\mapsto G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w})$ is
continuous for fixed $\boldsymbol{v}\in V_{2}$. Next, we verify that
$\boldsymbol{v}\mapsto G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w})$ is
continuous for fixed $\boldsymbol{w}\in V_{2}$. Let
$\boldsymbol{v},\boldsymbol{w}\in V_{2}$, and $(\boldsymbol{v_{k}})_{k}$ be a
sequence in $V_{2}$ with $\boldsymbol{v_{k}}\to\boldsymbol{v}\in V_{2}$. The
convergence of $\nabla\boldsymbol{v_{k}}\to\nabla\boldsymbol{v}$ in
$L^{2}(\Omega)^{N\times N}$ implies with the Riesz-Fischer theorem that we
find a subsequence $(\boldsymbol{v_{k}})_{k}$ with
$\nabla\boldsymbol{v_{k}}\to\nabla\boldsymbol{v}$ almost everywhere. Thus, we
have the pointwise convergence of the integrand in equation $(3.2)$.

We obtained an upper bound for $\langle
G_{1}^{\prime}(\boldsymbol{v_{k}})\boldsymbol{w},\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}$
with relation $(\ref{UpperG1})$ independently of $\boldsymbol{v_{k}}$.
Therefore, we apply the dominated convergence theorem to conclude

$\displaystyle G_{1}^{\prime}(\boldsymbol{v_{k}})\boldsymbol{w}\to
G_{1}^{\prime}(\boldsymbol{v})\boldsymbol{w}.$

We already proved that $\boldsymbol{w}\mapsto
G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w})$ is continuous for fixed
$\boldsymbol{v}\in V_{2}$. Thus, the mapping
$(\boldsymbol{v},\boldsymbol{w})\mapsto
G_{1}^{\prime}(\boldsymbol{v};\boldsymbol{w})$ is continuous. This implies the
Fréchet differentiability. ∎

###  3.2. Regularized p-Stokes equations

The solution of the $p$-Stokes equations is in $V_{p}$, see Theorem 2.8.
Differentiability is only possible in $V_{2}$, see Remark 3.2 and Theorem 3.3.
Hence, we need a smoother solution that is in $V_{2}$. Thus, we add a
regularization term to consider the problem on $V_{2}$. Let $G:V_{2}\to
V_{2}^{*}$,

$\displaystyle\langle
G(\boldsymbol{v}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}=\langle
A(\boldsymbol{v}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}+(\rho\boldsymbol{g},\boldsymbol{\phi})+\mu_{0}(\nabla\boldsymbol{v},\nabla\boldsymbol{\phi})\text{,
}\quad\boldsymbol{\phi}\in V_{2}$

with $\mu_{0}>0$. We call a function $\boldsymbol{v}\in V_{2}$ with

(12) $\displaystyle\langle
G(\boldsymbol{v}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}=0\quad\text{ for
all }\boldsymbol{\phi}\in V_{2}$

a weak solution of the regularized $p$-Stokes equations. We prove the
existence and uniqueness of solutions for the regularized $p$-Stokes
equations. For that purpose, we verify the continuity of $G:V_{2}\to
V_{2}^{*}$. In fact, we prove Lipschitz continuity.

######  Lemma 3.4.

The operator $G:V_{2}\to V_{2}^{*}$ is Lipschitz continuous.

###### Proof.

The mapping
$\boldsymbol{v}\mapsto\mu_{0}(\nabla\boldsymbol{v},\nabla\boldsymbol{\phi})$
is Lipschitz continuous for fixed $\boldsymbol{\phi}\in V_{2}$. In Lemma 2.5,
we proved that $A:V_{2}\to V_{2}^{*}$ is Lipschitz continuous. Thus, the
operator $G$ is Lipschitz continuous. ∎

Now the existence result follows easily:

######  Theorem 3.5.

There exists exactly one weak solution for the regularized $p$-Stokes
equations, see equation $(\ref{WeakRegulaizedSolution})$.

###### Proof.

We already concluded strict monotonicity in Theorem 2.8 in $V_{p}$. Thus, the
strict monotonicity follows for $A:V_{2}\to V_{2}$. In Lemma 2.7, we proved
coercivity in the $V_{p}$-norm. This leads to $\langle
A\boldsymbol{v},\boldsymbol{v}\rangle_{V_{2}^{*},V_{2}}\geq 0$. The strict
monotonicity and coercivity of $A+\mu_{0}$ follows with

$\displaystyle\mu_{0}(\nabla\boldsymbol{v},\nabla\boldsymbol{v})=\mu_{0}\|\nabla\boldsymbol{v}\|_{V_{2}}^{2}\text{,
}\quad\mu_{0}(\nabla(\boldsymbol{v}-\boldsymbol{w}),\nabla(\boldsymbol{v}-\boldsymbol{w}))=\mu_{0}\|\nabla(\boldsymbol{v}-\boldsymbol{w})\|_{V_{2}}^{2}>0$

for all $\boldsymbol{v},\boldsymbol{w}\in V_{2}$. Lemma 3.4 yields continuity.
We apply the Browder-Minty Theorem 2.2 to conclude that there exists exactly
one weak solution for $(\ref{WeakRegulaizedSolution})$. ∎

We apply the same arguments for the projection into divergence-free spaces as
in section 2.1. Hence, we obtain a unique pressure $\pi\in L^{2}_{0}(\Omega)$
for the regularized $p$-Stokes equations.

We consider the derivative of $G$ now:

######  Theorem 3.6.

Let $\boldsymbol{v},\boldsymbol{w}\in V_{2}$. The directional derivative has
the form

(13) $\displaystyle\begin{split}&\quad\langle
G^{\prime}(\boldsymbol{v})\boldsymbol{w},\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}\\\
&=\int_{\Omega}(p-2)B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-4)/2}(D\boldsymbol{v}:D\boldsymbol{w})\,(D\boldsymbol{v}:\nabla\boldsymbol{\phi})\,dx\\\
&\quad+\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}D\boldsymbol{w}:\nabla\boldsymbol{\phi}\,dx\\\
&\quad+\int_{\Gamma_{b}}(s-2)\tau\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(s-4)/2}(\boldsymbol{v}\cdot\boldsymbol{w})\,(\boldsymbol{v}\cdot\boldsymbol{\phi})\,ds\\\
&\quad+\int_{\Gamma_{b}}\tau\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(s-2)/2}\boldsymbol{w}\cdot\boldsymbol{\phi}\,ds+\mu_{0}\int_{\Omega}\nabla\boldsymbol{w}:\nabla\boldsymbol{\phi}\,dx\text{,
}\quad\boldsymbol{\phi}\in V_{2}.\end{split}$

Furthermore, the operator $G:V_{2}\to V_{2}^{*}$ is Fréchet differentiable.

###### Proof.

We calculated the directional derivatives for $G_{1}$ in Lemma 3.1. The
directional derivatives for the boundary term follow in the same way if we use
as the integration area $\Gamma_{b}$ instead of $\Omega$ and $V,\Phi$ with
$V_{ij}:=I_{ij}v_{j}$ and $\Phi_{ij}:=I_{ij}\phi_{j}$ instead of
$D\boldsymbol{v}$ and $D\boldsymbol{\phi}$. With this notation, we relate the
vector-valued expressions to the matrix-valued expressions. The arguments as
in Theorem 3.3 and the trace operator yield the Fréchet differentiability of
the boundary term. The last summand in equation $(\ref{Derivative})$ follows
trivially due to the linearity of
$\boldsymbol{v}\mapsto\mu_{0}(\nabla\boldsymbol{v},\nabla\boldsymbol{\phi})$,
$\boldsymbol{\phi}\in V_{2}$. ∎

##  4\. Infinite-dimensional Newton’s Method

In this section, we state Newton’s method in infinite dimensions as in [19]
and prove local convergence.

Newton’s method: Choose $\boldsymbol{v_{0}}\in V_{2}$ sufficiently close to
the solution $\boldsymbol{\overline{v}}$ of $G(\boldsymbol{\overline{v}})=0$.
For $k=0,1,2,...:$
Obtain $\boldsymbol{s_{k}}$ by solving

(14) $\displaystyle
G^{\prime}(\boldsymbol{v_{k}})\boldsymbol{s_{k}}=-G(\boldsymbol{v_{k}}),$

and set $\boldsymbol{v_{k+1}}:=\boldsymbol{v_{k}}+\boldsymbol{s_{k}}$.

######  Theorem 4.1.

Let $\bar{\boldsymbol{v}}$ satisfy $G(\bar{\boldsymbol{v}})=0$. Let the start
point $\boldsymbol{v_{0}}$ be sufficiently close to the solution
$\bar{\boldsymbol{v}}$ and $(\boldsymbol{v_{k}})_{k}$ be generated by the
infinite-dimensional Newton’s method, see equation $(\ref{NewtonAlg})$. Set
$\boldsymbol{d_{k}}:=\boldsymbol{v_{k}}-\bar{\boldsymbol{v}}$. Let further the
_regularity condition_

(15)
$\displaystyle\|G^{\prime}(\boldsymbol{v_{k}})^{-1}\|_{\mathcal{L}(V_{2},V_{2}^{*})}\leq
C$

independently of $k$ and the _approximation condition_

(16)
$\displaystyle\|G(\bar{\boldsymbol{v}}+\boldsymbol{d_{k}})-G(\bar{\boldsymbol{v}})-G^{\prime}(\boldsymbol{v_{k}})\boldsymbol{d_{k}}\|_{V_{2}^{*}}=o(\|\boldsymbol{d_{k}}\|_{V_{2}})$

for $\|\boldsymbol{d_{k}}\|_{V_{2}}\to 0$ be fulfilled for all
$k\in\mathbb{N}$. Then $\boldsymbol{v_{k}}$ converges $q$-superlinearly to
$\bar{\boldsymbol{v}}$.

The problem to calculate a Newton step, see equation $(\ref{NewtonAlg})$, is
linear because we handle the linear problem
$(\boldsymbol{w},\boldsymbol{\phi})\mapsto\langle
G^{\prime}(\boldsymbol{v})\boldsymbol{w},\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}$
for all $\boldsymbol{w},\boldsymbol{\phi}\in V_{2}$. Hence, we can apply Lax-
Milgram’s Lemma. We need to verify the continuity and coercivity of this
linear problem to obtain the existence and uniqueness of a solution to that
problem. Additionally, we prove the conditions in Theorem 4.1 to apply the
infinite-dimensional Newton’s method. The function $G^{\prime}$ is continuous.
For continuously Fréchet differentiable functions the approximation condition
$(\ref{Convergence})$ is fulfilled, see [19, section 2.4.3]. Hence, it remains
to show the regularity condition.

######  Lemma 4.2 (Regularity condition).

Let $\boldsymbol{v}\in V_{2}$. There exists exactly one solution
$\boldsymbol{w}=\boldsymbol{w}(\boldsymbol{v})$ such that equation
$(\ref{NewtonAlg})$ is fullfiled. Moreover, the regularity condition, see
equation $(\ref{RegCond})$, is true with $C:=1/\mu_{0}$.

###### Proof.

Let $\boldsymbol{v}\in V_{2}$. We show the continuity of
$(\boldsymbol{w},\boldsymbol{\phi})\mapsto\langle
G^{\prime}(\boldsymbol{v})\boldsymbol{w},\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}$
for all $\boldsymbol{w},\boldsymbol{\phi}\in V_{2}$. Since the arguments in
the calculation of relation $(\ref{Bound})$ are also valid for the boundary
term, we find $\tilde{C}\in\mathbb{R}$ with the trace operator, see [19,
Theorem 1.12], such that

$\displaystyle|\langle
G^{\prime}(\boldsymbol{v})\boldsymbol{w},\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}|\leq\tilde{C}(\delta^{s-2}\|\tau\|_{L^{\infty}(\Gamma_{b})}+\delta^{p-2}\|B\|_{L^{\infty}(\Omega)}+\mu_{0})\|\boldsymbol{w}\|_{V_{2}}\|\boldsymbol{\phi}\|_{V_{2}}.$

We show the coercivity next. Let $\boldsymbol{w}\in V_{2}$. The last summand
of the directional derivative of $G^{\prime}$, see equation
$(\ref{Derivative})$, is just

$\displaystyle\mu_{0}(\nabla\boldsymbol{w},\nabla\boldsymbol{w})=\mu_{0}\|\boldsymbol{w}\|_{V_{2}}.$

We use $D\boldsymbol{w}:\nabla\boldsymbol{v}=D\boldsymbol{w}:D\boldsymbol{v}$
for arbitrary $\boldsymbol{v},\boldsymbol{w}\in V_{2}$, see equation
$(\ref{Transpose})$, to conclude

$\displaystyle\quad\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}D\boldsymbol{w}:\nabla\boldsymbol{w}\,dx$
$\displaystyle\quad+\int_{\Omega}(p-2)B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-4)/2}(D\boldsymbol{v}:D\boldsymbol{w})\,(D\boldsymbol{v}:\nabla\boldsymbol{w)}\,dx$
$\displaystyle=\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\bigg{(}|D\boldsymbol{w}|^{2}+(p-2)\frac{(D\boldsymbol{v}:D\boldsymbol{w})^{2}}{|D\boldsymbol{v}|^{2}+\delta^{2}}\bigg{)}\,dx$
$\displaystyle\geq\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\bigg{(}|D\boldsymbol{w}|^{2}+(p-2)\frac{|D\boldsymbol{v}|^{2}\,|D\boldsymbol{w}|^{2}}{|D\boldsymbol{v}|^{2}}\bigg{)}\,dx$
$\displaystyle=\int_{\Omega}B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}(p-1)|D\boldsymbol{w}|^{2}\,dx\geq
0.$

The same arguments for the boundary term yield

$\displaystyle\int_{\Gamma_{b}}\tau\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(s-2)/2}|\boldsymbol{w}|^{2}+(s-2)\tau\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(s-4)/2}(\boldsymbol{v}\cdot\boldsymbol{w})\,(\boldsymbol{v}\cdot\boldsymbol{w})\,ds\geq
0.$

Hence, the conditions for applying Lax-Milgram’s Lemma are validated, and a
unique solution for each Newton step, see equation $(\ref{NewtonAlg})$,
exists. The coercivity constant is $\mu_{0}$. This implies the uniform bound

$\displaystyle\|(G^{\prime}(\boldsymbol{v}))^{-1}\|_{\mathcal{L}(V_{2},V_{2}^{*})}\leq
1/\mu_{0}.$

∎

Because we proved the well-posedness, the regularity condition, and the
approximation condition, we conclude:

######  Corollary 4.3 (Local $q$-superlinear convergence).

The infinite-dimensional Newton’s method converges locally $q$-superlinearly
to a solution for the regularized $p$-Stokes equations in the divergence-free
space, see equation $(\ref{WeakRegulaizedSolution})$.

##  5\. Global q-superlinear convergence

In this section, we globalize Newton’s method with a step size control and
prove that this modification merges into the local $q$-superlinear convergent
Newton’s method, see equation $(\ref{NewtonAlg})$. To globalize Newton’s
method, we formulate a minimization problem using the fact that solving the
$p$-Stokes equations is equivalent to minimizing a convex functional.

######  Remark 5.1.

The intuitive approach may be to determine step sizes by minimizing
$\|G(\boldsymbol{v})\|_{V_{2}^{*}}^{2}$. This leads us either to a finite-
dimensional consideration with the Euclidean norm or to calculating the
residual norm, which corresponds in solving an additional Stokes problem for
each iteration. Because we stick to the optimize then discretize approach, we
use an alternative described in the following.

###  5.1. The convex functional

In this subsection, we introduce the convex functional and verify that the
directional derivative is the variational formulation tested with the
directional derivative. In [4], the convex functional was introduced for
$\mu_{0}=0$ and Dirichlet boundary conditions. In [2], it was introduced for
$\mu_{0}=0$ and $\delta=0$ for more general boundary conditions. We use those
formulations for our convex functional with $\mu_{0},\delta>0$ or
$(\mu_{0},\delta)=(0,0)$.

######  Lemma 5.2.

Let $\mu_{0},\delta\in(0,\infty)$ or $(\mu_{0},\delta)$=0 and
$J_{\mu_{0},\delta}:V_{2}\to\mathbb{R}$,

$\displaystyle\begin{split}J_{\mu_{0},\delta}(\boldsymbol{v})=&\int_{\Omega}\frac{B}{p}\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{p/2}\,dx+\int_{\Gamma_{b}}\frac{\tau}{s}\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{s/2}\,ds+\mu_{0}(\nabla\boldsymbol{v},\nabla\boldsymbol{v})+(\rho\boldsymbol{g},\boldsymbol{v}).\end{split}$

The functional $J_{\mu_{0},\delta}$ is Fréchet differentiable with

$\displaystyle
J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v})\boldsymbol{w}=\langle
A(\boldsymbol{v}),\boldsymbol{w}\rangle_{V_{2}^{*},V_{2}}+\mu_{0}(\nabla\boldsymbol{v},\nabla\boldsymbol{w})+(\rho\boldsymbol{g},\boldsymbol{w}).$

###### Proof.

The derivative of the third summand of $J_{\mu_{0},\delta}$ is trivial. Thus,
we consider $\mu_{0}=0$. Let $\boldsymbol{v},\boldsymbol{w}\in V_{2}$,
$t\in(0,1)$. We insert the definition to obtain

$\displaystyle\quad\frac{1}{t}\big{(}J_{0,\delta}(\boldsymbol{v}+t\boldsymbol{w})-J_{0,\delta}(\boldsymbol{v})\big{)}$
$\displaystyle=\int_{\Omega}\frac{B}{pt}\Big{(}\big{(}|D(\boldsymbol{v}+t\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{p/2}-\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{p/2}\Big{)}\,dx$
$\displaystyle\quad+\int_{\Gamma_{b}}\frac{\tau}{st}\Big{(}\big{(}|\boldsymbol{v}+t\boldsymbol{w}|^{2}+\delta^{2}\big{)}^{p/2}-\big{(}|\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{p/2}\Big{)}\,ds+(\rho\boldsymbol{g},\boldsymbol{w}).$

We consider the first summand and define the auxiliary function
$f:\Omega\times\mathbb{R}\to\mathbb{R}$,

$\displaystyle
f(x,t)=\big{(}|D(\boldsymbol{v}(x)+t\boldsymbol{w}(x))|^{2}+\delta^{2}\big{)}^{p/2}.$

This definition yields for the derivative

$\displaystyle\frac{\partial f}{\partial t}(\cdot,t)=$ $\displaystyle
p\big{(}|D(\boldsymbol{v}+t\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{(p-2)/2}(D\boldsymbol{v}:D\boldsymbol{w}+t|D\boldsymbol{w}|^{2}).$

We apply the Taylor expansion to obtain with $\xi:\Omega\to[0,t]$,

$\displaystyle\quad\int_{\Omega}\frac{B}{pt}\Big{(}\big{(}|D(\boldsymbol{v}+t\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{p/2}-\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{p/2}\Big{)}\,dx$
$\displaystyle=\int_{\Omega}\frac{B}{pt}\big{(}f(x,0)+f^{\prime}(x,\xi(x))t-f(x,0)\big{)}\,dx$
$\displaystyle=\int_{\Omega}\frac{B}{pt}p\big{(}|D(\boldsymbol{v}+\xi\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{(p-2)/2}(D\boldsymbol{v}:D\boldsymbol{w}+\xi|D\boldsymbol{w}|^{2})t\,dx$
$\displaystyle=\int_{\Omega}B\big{(}|D(\boldsymbol{v}+\xi\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{(p-2)/2}(D\boldsymbol{v}:D\boldsymbol{w}+\xi|D\boldsymbol{w}|^{2})\,dx.$

Due to $\xi\leq t$, the pointwise convergence of the integrand for $t\to 0$
can be easily seen:

$\displaystyle
B\big{(}|D(\boldsymbol{v}+\xi\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{(p-2)/2}(D\boldsymbol{v}:D\boldsymbol{w}+\xi|D\boldsymbol{w}|^{2})$
$\displaystyle\to
B\big{(}|D\boldsymbol{v}|^{2}+\delta^{2}\big{)}^{(p-2)/2}D\boldsymbol{v}:D\boldsymbol{w}$
$\displaystyle=BS^{p}(D\boldsymbol{v}):D\boldsymbol{w}.$

Furthermore, we obtain the upper bound

$\displaystyle\quad
B\big{(}|D(\boldsymbol{v}+\xi\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{(p-2)/2}(D\boldsymbol{v}:D\boldsymbol{w}+\xi|D\boldsymbol{w}|^{2})$
$\displaystyle=B(|D(\boldsymbol{v}+\xi\boldsymbol{w})|^{2}+\delta^{2})^{(p-2)/2}(D(\boldsymbol{v}+\xi\boldsymbol{w}):D\boldsymbol{w})$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}(|D(\boldsymbol{v}+\xi\boldsymbol{w})|)^{p-1}|D\boldsymbol{w}|$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}(|D\boldsymbol{v}|+|D\boldsymbol{w}|)^{p-1}|D\boldsymbol{w}|.$

Thereby, we can use dominated convergence independently of $\delta\geq 0$ to
conclude

$\displaystyle\lim_{t\to 0}$
$\displaystyle\int_{\Omega}B\big{(}|D(\boldsymbol{v}+\xi\boldsymbol{w})|^{2}+\delta^{2}\big{)}^{(p-2)/2}(D\boldsymbol{v}:D\boldsymbol{w}+\xi|D\boldsymbol{w}|^{2})\,dx=\int_{\Omega}BS^{p}(D\boldsymbol{v}):D\boldsymbol{w}\,dx.$

The arguments are identical for the second summand. Thus, we calculated the
directional derivative of $J_{0,\delta}$. The linearity with respect to the
direction $\boldsymbol{w}$ implies the Gâteaux differentiability. The
continuity of $A$, see Lemma 2.5, and the continuity of the second summand
imply the Fréchet differentiability:

$\displaystyle J_{0,\delta}^{\prime}(\boldsymbol{v})\boldsymbol{w}=\langle
A(\boldsymbol{v}),\boldsymbol{w}\rangle_{V_{2}^{*},V_{2}}+(\nabla\boldsymbol{v},\nabla\boldsymbol{w})+(\rho\boldsymbol{g},\boldsymbol{w})\in
V_{p}.$

∎

We relate the $p$-Stokes equations with the convex functional:

######  Lemma 5.3.

We have for the case $\mu_{0},\delta\in(0,1)$ and the case $\mu_{0}=0=\delta$

$\displaystyle J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v})=0\Leftrightarrow
J_{\mu_{0},\delta}(\boldsymbol{v})=\min_{\boldsymbol{u}\in
V_{p}}J_{\mu_{0},\delta}(\boldsymbol{u}).$

###### Proof.

We verified in Lemma 15 positive definitness of $J_{\mu_{0},\delta}$ for
$(\mu_{0},\delta)\in(0,\infty)^{2}$. This implies strict convexity. The strict
convexity for $(\mu_{0},\delta)=(0,0)$ was proved in [2]. Because the
necessary first-order optimality condition is sufficient for strict convex
functions, the claim follows. ∎

######  Remark 5.4.

We only verified existence and uniqueness of a solution for
$\mu_{0},\delta>0$. In [2], the case $\mu_{0}=0=\delta$ was considered. Thus,
we exclude $\mu_{0}=0$ and $\delta>0$, and $\mu_{0}>0$ and $\delta=0$ to
shorten the arguments. Moreover, these cases are not necessary for the
upcoming analysis.

###  5.2. Global q-superlinear convergence

In this subsection, we verify global convergence. In a second step, we verify
$q$-superlinar convergence. We obtain global convergence by employing a step
size control. This step size control considers the problem to minimize
$J_{\mu_{0},\delta}$, see Lemma $(\ref{ConvexFunction})$. Conditions for a
globally convergent method are given in [19, Theorem 5.2]. We verify these
conditions for given $\boldsymbol{v_{k}}\in V_{2}$ and search direction
$0\neq\boldsymbol{w_{k}}\in V_{2}$:

  1. a)

We have a descent direction:
$J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\boldsymbol{w_{k}}<0$.

  2. b)

The angle condition is fulfilled:

$\displaystyle\eta\|J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\|_{V_{2}^{*}}\|\boldsymbol{w_{k}}\|_{V_{2}}\leq-
J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\boldsymbol{w_{k}}\quad\text{
independently of }k.$

  3. c)

The step sizes are admissible:

$\displaystyle
J_{\mu_{0},\delta}(\boldsymbol{v_{k}}+\alpha_{k}\boldsymbol{w_{k}})$
$\displaystyle<J_{\mu_{0},\delta}(\boldsymbol{v_{k}})\quad\text{ for all
}k\text{ and }$ $\displaystyle
J_{\mu_{0},\delta}(\boldsymbol{v_{k}}+\alpha_{k}+\boldsymbol{w_{k}})$
$\displaystyle-J_{\mu_{0},\delta}(\boldsymbol{v_{k}})\to 0\quad\text{ for
}k\to\infty$ $\displaystyle\text{ imply
}\frac{J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\boldsymbol{w_{k}}}{\|\boldsymbol{w_{k}}\|_{V_{2}}}\to
0\quad\text{ for }k\to\infty.$

The first two statements are easy to verify:

######  Lemma 5.5.

Newton steps fulfill the angle condition for the functional
$J_{\mu_{0},\delta}$.

###### Proof.

We calculated in Lemma 4.2 the coercivity constant $\mu_{0}$ and the
continuity constant $C$ independent of $k$ such that

(17)
$\displaystyle\mu_{0}\leq\|J_{\mu_{0},\delta}^{\prime\prime}(\boldsymbol{v_{k}})\|_{\mathcal{L}(V_{2},V_{2}^{*})}\leq
C.$

In [21, Theorem 1.2], the proof of the angle condition is done with
$\eta\geq\mu_{0}/C$ for the finite-dimensional case. The proof for the
infinite-dimensional case is identical. ∎

We obtain admissible step sizes:

######  Lemma 5.6.

Let $\beta\in(0,1)$, $\gamma\in(0,1/2)$, $\boldsymbol{v_{k}}\in V_{2}$ the
point, and $\boldsymbol{s_{k}}\in V_{2}$ the direction. Determine the Armijo
step sizes, see [19, section 2.2.1.1], by finding the biggest
$\alpha_{k}\in\\{1,\beta,\beta^{2},...\\}$ such that

(18) $\displaystyle
J_{\mu_{0},\delta}(\boldsymbol{v_{k}}+\alpha_{k}\boldsymbol{s_{k}})-J_{\mu_{0},\delta}(\boldsymbol{v_{k}})\leq\alpha_{k}\gamma
J^{\prime}_{\mu_{0},\delta}(\boldsymbol{v_{k}})\boldsymbol{s_{k}}.$

For Newton’s method, these Armijo step sizes are admissible.

###### Proof.

We want to apply [19, Lemma 2.3], which states that $J_{\mu_{0},\delta}$ has
to be uniformly continuous, and the step sizes have to be large enough. We
verify these conditions. The function $J^{\prime}_{\mu_{0},\delta}$ is
Lipschitz-continuous, see Lemma 3.4, and therefore uniformly continuous. We
already proved that the directions $\boldsymbol{w_{k}}$ are descent
directions. It remains to show that the descent directions are not too short.
We know with the Newton steps and the right inequality in relation
$(\ref{BoundJprimeprime})$

$\displaystyle\|J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\|_{V_{2}^{*}}=\|J_{\mu_{0},\delta}^{\prime\prime}(\boldsymbol{v_{k}})\boldsymbol{w_{k}}\|_{V_{2}^{*}}\leq
C\|\boldsymbol{w_{k}}\|_{V_{2}}.$

Thus, we obtain

$\displaystyle\|\boldsymbol{w_{k}}\|_{V_{2}}\geq\frac{1}{C}\|J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\|_{V_{2}^{*}}\geq\frac{1}{C}\frac{-J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{k}})\boldsymbol{w_{k}}}{\|\boldsymbol{w_{k}}\|_{V_{2}}}.$

Therefore, the step size is bounded from below by the monotonically increasing
function $\phi:[0,\infty)\to[0,\infty)$, $\phi(t)=t/C$ with $\phi(t)>0$ for
$t>0$. Hence, the step sizes are large enough and we can apply [19, Lemma
2.3]. ∎

######  Corollary 5.7.

Newton’s method with Armijo step sizes is globally convergent.

###### Proof.

We already explained why the step sizes are admissible, see Lemma 5.6.
Moreover, the directions $\boldsymbol{w_{k}}$ are descent directions.
Furthermore, we validated the angle condition. The angle condition implies the
admissibility of search directions, see [19, Lemma 2.1]. Trivially, the
functional $J_{\mu_{0},\delta}$ is bounded from below. Thus, all conditions
for [19, Theorem 2.2] are fulfilled, and the claim follows. ∎

######  Remark 5.8.

The $q$-superlinear convergence follows if we verify that the Armijo rule
chooses step sizes of length $1$ after some iterations. The finite-dimensional
case for quasi-Newton methods is given in [22, Proposition 1.15]. The
infinite-dimensional proof with the Hessian is identical because
$J^{\prime\prime}_{\mu_{0},\delta}$ is positive definite and symmetric.

##  6\. Convergence to the non regularized problem

In this section, we argue that the solution
$\boldsymbol{v}=\boldsymbol{v_{\mu_{0},\delta}}$ converges to
$\boldsymbol{v_{0,0}}$ in $L^{p}(\Omega)^{N}$ for vanishing $\mu_{0}$ and
$\delta$. This is important because we have $\delta=0$ and $\mu_{0}=0$ for
glaciological applications. We do this with the following intermediate steps:

  * a1)

We show that the functional vanishes in $V_{q}^{*}$:

$\displaystyle\|G_{0,0}\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{q}^{*}}\leq\epsilon\quad\text{
for all }\mu_{0}\in(0,\hat{\mu})\text{ and }\delta\in(0,\hat{\delta}).$

  * a2)

We bound the solution in $V_{p}$:

$\displaystyle\|\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{p}}\leq c\quad\text{ for
all }\mu_{0},\delta\in(0,1),$

  * b)

We show that the functional vanishes in $V_{p}^{*}$:

$\displaystyle\|G_{0,0}\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{p}^{*}}\leq\epsilon\quad\text{
for all }\mu_{0}\in(0,\hat{\mu})\text{ and }\delta\in(0,\hat{\delta}).$

Furthermore, we prove
$\boldsymbol{v_{\mu_{0},\delta}}\to\boldsymbol{v_{0,0}}\in V_{2}$, if
$\boldsymbol{v_{0,0}}\in V_{2}$ is fulfilled.

We do not use the formulation $(\mu_{0},\delta)\downarrow(0,0)$ because we did
not prove existence and uniqueness of a solution for $\mu_{0}>0$ and
$\delta=0$. Thus, the existence and uniqueness of $\boldsymbol{v_{\mu_{0},0}}$
is not clear for $\mu_{0}>0$.

######  Lemma 6.1.

Let $\boldsymbol{v_{\mu_{0},\delta}}$ be the solution of
$G_{\mu_{0},\delta}(\boldsymbol{v}):=G(\boldsymbol{v})=0$ for
$\mu_{0},\delta\in(0,1)$ with the variables $\mu_{0}$ and $\delta$ as in the
definition of the operator $G$. For all $\epsilon>0$, there exist
$\hat{\mu},\hat{\delta}\in(0,1)$ with

$\displaystyle\|G_{0,0}\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{q}^{*}}\leq\epsilon$

for all $\mu_{0}\in(0,\hat{\mu})$ and $\delta\in(0,\hat{\delta})$. Moreover,
there exists $c\in\mathbb{R}$ with

$\displaystyle\|\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{p}}\leq c\quad\text{ for
all }\mu_{0},\delta\in(0,1).$

###### Proof.

We set $\boldsymbol{v_{k}}:=\boldsymbol{v_{\mu_{0},\delta}}$. Let
$\boldsymbol{\phi}\in V_{q}$. We conclude

$\displaystyle\quad|\langle
G_{0,0}(\boldsymbol{v_{k}}),\boldsymbol{\phi}\rangle_{V_{p}^{*},V_{q}}|$
$\displaystyle=|\langle
G_{\mu_{0},\delta}(\boldsymbol{v_{k}}),\boldsymbol{\phi}\rangle_{V_{p}^{*},V_{q}}-\mu_{0}(\nabla\boldsymbol{v_{k}},\nabla\boldsymbol{\phi})$
$\displaystyle\quad+\int_{\Omega}B\Big{(}|D\boldsymbol{v_{k}}|^{p-2}-\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\Big{)}D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx$
$\displaystyle\quad+\int_{\Gamma_{b}}\tau\Big{(}|\boldsymbol{v_{k}}|^{p-2}-\big{(}|\boldsymbol{v_{k}}|+\delta^{2}\big{)}^{(p-2)/2}\Big{)}\boldsymbol{v_{k}}\cdot\boldsymbol{\phi}\,ds\bigg{|}$
(19)
$\displaystyle\leq\mu_{0}\|\boldsymbol{v_{k}}\|_{V_{p}}\|\boldsymbol{\phi}\|_{V_{q}}+\bigg{|}\int_{\Omega}B\Big{(}|D\boldsymbol{v_{k}}|^{p-2}-\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\Big{)}D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
$\displaystyle\quad+\bigg{|}\int_{\Gamma_{b}}\tau\Big{(}|\boldsymbol{v_{k}}|^{p-2}-\big{(}|\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\Big{)}\boldsymbol{v_{k}}\cdot\boldsymbol{\phi}\,ds\bigg{|}.$

First, we consider the second summand. We split the integral into two
components and verify that both can be made arbitrarily small for small enough
$\delta$ and $\mu_{0}$. Let $\epsilon>0$. Set
$M_{\epsilon}:=\\{\boldsymbol{x}\in\Omega\text{;
}|D\boldsymbol{v_{k}}(\boldsymbol{x})|<\epsilon\\}$. We obtain

$\displaystyle\quad\bigg{|}\int_{M_{\epsilon}}B\Big{(}|D\boldsymbol{v_{k}}|^{p-2}-\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\Big{)}D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}\int_{M_{\epsilon}}\Big{(}|D\boldsymbol{v_{k}}|^{p-2}+\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\Big{)}|D\boldsymbol{v_{k}}|\,|\nabla\boldsymbol{\phi}|\,dx$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}\int_{M_{\epsilon}}2|D\boldsymbol{v_{k}}|^{p-1}|\nabla\boldsymbol{\phi}|\,dx$
$\displaystyle\leq
2\epsilon^{p-1}\|B\|_{L^{\infty}(\Omega)}\|\boldsymbol{\phi}\|_{V_{q}}\|1\|_{V_{p}}.$

We consider the complement of $M_{\epsilon}$ next. Let
$\boldsymbol{x}\in\Omega\setminus M_{\epsilon}$. We define
$f:\Omega\times\mathbb{R}\to\mathbb{R}$,

$\displaystyle
f(x,t)=\big{(}|D\boldsymbol{v_{k}}(x)|^{2}+t^{2}\big{)}^{(p-2)/2}\text{,
}\quad
f^{\prime}(x,t)=(p-2)\big{(}|D\boldsymbol{v_{k}}(x)|^{2}+t^{2}\big{)}^{(p-4)/2}t.$

A Taylor expansion yields with $\xi:\Omega\setminus M_{\epsilon}\to[0,\delta]$

$\displaystyle\quad\bigg{|}\int_{\Omega\setminus
M_{\epsilon}}B\Big{(}|D\boldsymbol{v_{k}}|^{p-2}-\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\Big{)}D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
$\displaystyle=\bigg{|}\int_{\Omega\setminus
M_{\epsilon}}B(f(\cdot,0)-f(\cdot,0)-f^{\prime}(\cdot,\xi)\delta)D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
$\displaystyle=\bigg{|}\int_{\Omega\setminus
M_{\epsilon}}B(2-p)\big{(}|D\boldsymbol{v_{k}}|^{2}+\xi^{2}\big{)}^{(p-4)/2}\xi\delta
D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
$\displaystyle\leq(2-p)\|B\|_{L^{\infty}(\Omega)}\int_{\Omega}\epsilon^{p-4}\delta^{2}|D\boldsymbol{v_{k}}|\,|\nabla\boldsymbol{\phi}|\,dx$
$\displaystyle\leq(2-p)\|B\|_{L^{\infty}(\Omega)}\epsilon^{p-4}\delta^{2}\|\boldsymbol{v_{k}}\|_{V_{p}}\|\boldsymbol{\phi}\|_{V_{q}}.$

Thus, we conclude

$\displaystyle\quad\bigg{|}\int_{\Omega}B\bigg{(}|D\boldsymbol{v_{k}}|^{p-2}-\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}\bigg{)}D\boldsymbol{v_{k}}:\nabla\boldsymbol{\phi}\,dx\bigg{|}$
$\displaystyle\leq
2\epsilon^{p-1}\|B\|_{L^{\infty}(\Omega)}\|\boldsymbol{\phi}\|_{V_{q}}\|1\|_{V_{p}}+(2-p)\|B\|_{L^{\infty}(\Omega)}\epsilon^{p-4}\delta^{2}\|\boldsymbol{v_{k}}\|_{V_{p}}\|\boldsymbol{\phi}\|_{V_{q}}.$

For small enough $\epsilon$, the first summand becomes arbitrarily small. For
bounded $\|\boldsymbol{v_{k}}\|_{V_{p}}$, we choose $\delta$ small enough to
have an arbitrarily small second summand. It remains to find an upper bound
for $\|\boldsymbol{v_{k}}\|_{V_{p}}$. The function $\boldsymbol{v_{k}}$ is the
solution of

$\displaystyle\langle
G_{k}(\boldsymbol{v_{k}}),\boldsymbol{\phi}\rangle_{V_{2}^{*},V_{2}}=-(\rho\boldsymbol{g},\boldsymbol{\phi}).$

We set $\boldsymbol{\phi}:=\boldsymbol{v_{k}}$. Thus, we conclude with the
coercivity estimate $(\ref{CoerzEstimate})$ with $\tilde{c}\in\mathbb{R}$
independent of $\delta$ and $\mu_{0}$:

$\displaystyle\tilde{c}\|\boldsymbol{v_{k}}\|_{V_{p}}^{p}\leq\int_{\Omega}B\big{(}|D\boldsymbol{v_{k}}|^{2}+\delta^{2}\big{)}^{(p-2)/2}|D\boldsymbol{v_{k}}|^{2}\,dx\leq
c\|\rho\boldsymbol{g}\|_{L^{q}(\Omega)^{N}}\|\boldsymbol{v_{k}}\|_{V_{p}}.$

Thus, we obtain an upper bound for $\|\boldsymbol{v_{k}}\|_{V_{p}}$
independent of $k$.

The last summand on the right-hand side of inequality $(\ref{LastSummand})$ is
similar: Instead of $D\boldsymbol{v_{k}}$, $\nabla\boldsymbol{\phi}$, and $B$,
we have $\boldsymbol{v_{k}}$, $\boldsymbol{\phi}$, and $\tau$. We define the
set $\tilde{M}_{\epsilon}:=\\{\mathbf{x}\in\Gamma_{b}\text{;
}|\boldsymbol{v_{k}}(\boldsymbol{x})<\epsilon\\}$. The upcoming estimations
are the same with three exceptions: On $\tilde{M}_{\epsilon}$, we have an
upper bound with $\epsilon^{s-1}$ instead of $\epsilon^{p-1}$. On
$\Omega\setminus\tilde{M}_{\epsilon}$, we have $\epsilon^{s-4}$ instead of
$\epsilon^{p-4}$. Lastly, we apply the trace operator to obtain bounds on
$\Omega$ instead of $\Gamma_{b}$.

Hence, we have $G_{0,0}\boldsymbol{v_{k}}\to 0\in V_{q}^{*}$. ∎

######  Lemma 6.2.

For all $\epsilon>0$, there exist $\hat{\mu}$ and $\hat{\delta}$ with

$\displaystyle\|G_{0,0}\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{p}^{*}}\leq\epsilon\text{
for all }\mu_{0}\in(0,\hat{\mu})\text{ and }\delta\in(0,\hat{\delta}).$

###### Proof.

Let $(\mu_{0_{k}},\delta_{k})\in(0,1)^{\mathbb{N}}$ with
$(\mu_{0_{k}},\delta_{k})\downarrow(0,0)$. We set
$\boldsymbol{v_{k}}:=\boldsymbol{v_{\mu_{0_{k}},\delta_{k}}}$ and define the
linear operators $T_{k}:V_{p}\to\mathbb{R}$,

$\displaystyle T_{k}(\boldsymbol{\phi})=\langle
G_{0,0}\boldsymbol{v_{k}},\boldsymbol{\phi}\rangle_{V_{p}^{*},V_{p}}.$

We set $M:=\\{\boldsymbol{\phi}\in V_{p}\text{; }T_{k}(\boldsymbol{\phi})\to
0\text{ for }k\to\infty\\}$. Lemma 6.1 yields $V_{q}\subseteq M$. We will
extend this result to $M=V_{p}$. The set $\\{\boldsymbol{v}\in
C^{\infty}(\Omega)^{N}\text{; }\boldsymbol{v}|_{\Gamma_{d}=0}\\}$ is dense in
$V_{p}$. We know $\\{\boldsymbol{v}\in C^{\infty}(\Omega)^{N}\text{;
}\boldsymbol{v}|_{\Gamma_{d}}=0\\}\subseteq V_{q}\subseteq V_{p}$. Thus,
$V_{q}$ is a dense subset of $V_{p}$. Assume that $(T_{k})_{k}$ is uniformly
bounded, then $M$ is closed, see [23, Lemma 9.22]. This implies $M=V_{p}$ and
$G_{0,0}\boldsymbol{v_{k}}\to 0\in V_{p}^{*}$. Thus, it remains to verify that
$(T_{k})_{k}$ is uniformly bounded. We conclude for $\boldsymbol{\phi}\in
V_{p}$ with the dual exponent $q=p/(p-1)$, and $r:=p/(s-1)$, the dual exponent
$r^{\prime}$, and Hölder’s inequality

$\displaystyle|\langle T_{k},\boldsymbol{\phi}\rangle_{V_{p}^{*},V_{p}}|$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}\int_{\Omega}|D\boldsymbol{v_{k}}|^{p-1}|\nabla\boldsymbol{\phi}|\,dx+\|\tau\|_{L^{\infty}(\Gamma_{b})}\int_{\Gamma_{b}}|\boldsymbol{v_{k}}|^{s-1}|\boldsymbol{\phi}|\,ds$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}\bigg{(}\int_{\Omega}|D\boldsymbol{v_{k}}|^{p}\,dx\bigg{)}^{1/q}\|\boldsymbol{\phi}\|_{V_{p}}$
$\displaystyle\quad+\|\tau\|_{L^{\infty}(\Gamma_{b})}\bigg{(}\int_{\Gamma_{b}}|\boldsymbol{v_{k}}|^{p}\,ds\bigg{)}^{1/r}\|\boldsymbol{\phi}\|_{V_{r^{\prime}}}$
$\displaystyle\leq\|B\|_{L^{\infty}(\Omega)}\|\boldsymbol{v_{k}}\|_{V_{p}}^{p/q}\|\boldsymbol{\phi}\|_{V_{p}}+c\|\tau\|_{L^{\infty}(\Gamma_{b})}\|\boldsymbol{v_{k}}\|_{V_{p}}^{p/r}\|\boldsymbol{\phi}\|_{V_{p}}.$

The last inequality followed with $r>q$, which implies
$r^{\prime}<q^{\prime}=p$, and the trace operator. This is a uniform bound,
because $\|\boldsymbol{v_{k}}\|_{V_{p}}$ is uniformly bounded, see Lemma 6.1.
Thus, the claim follows. ∎

######  Corollary 6.3.

Let $(\mu_{0_{k}},\delta_{k})\in(0,1)^{\mathbb{N}}$ with
$(\mu_{0_{k}},\delta_{k})\downarrow(0,0)$. Set
$\boldsymbol{v_{k}}:=\boldsymbol{v_{\mu_{0_{k}},\delta_{k}}}$. It follows
$\boldsymbol{v_{k}}\to\boldsymbol{v_{0,0}}$ in $L^{p}(\Omega)^{N}$.

###### Proof.

The sequence $(\boldsymbol{v_{k}})_{k}$ is bounded, due to Lemma 6.1.
Consequently, there exists a subsequence $(\boldsymbol{v_{k_{m}}})_{m}$ with
$\boldsymbol{v_{k_{m}}}\rightharpoonup\boldsymbol{\overline{v}}\in V_{p}$.
Moreover, we know

$\displaystyle
A\boldsymbol{v_{k_{m}}}=G_{0,0}\boldsymbol{v_{k_{m}}}-\rho\boldsymbol{g}\to-\rho\boldsymbol{g}\in
V_{p}^{*},$

see Lemma 6.2. Thus, the decisive monotonicity trick, see [24, page 474],
already yields $A\boldsymbol{\overline{v}}=-\rho\boldsymbol{g}$, because $A$
is a continuous and monotone operator. The uniqueness of the solution for the
operator $A$ implies $\boldsymbol{v_{0,0}}=\boldsymbol{\overline{v}}$. We
conclude $\boldsymbol{v_{k_{m}}}\rightharpoonup\boldsymbol{v_{0,0}}\in V_{p}$.
This implies $\boldsymbol{v_{k_{m}}}\to\boldsymbol{v_{0,0}}$ in
$L^{p}(\Omega)^{N}$.

Assume: $\boldsymbol{v_{k}}\not\to\boldsymbol{v_{0,0}}$ in
$L^{p}(\Omega)^{N}$. Then, there exist $\epsilon>0$ and a subsequence
$(\boldsymbol{\tilde{v}_{k}})_{k}$ with

$\displaystyle\|\boldsymbol{\tilde{v}_{k}}-\boldsymbol{v_{0,0}}\|_{L^{p}(\Omega)^{N}}>\epsilon.$

But the boundedness of $(\boldsymbol{\tilde{v}_{k}})_{k}$ implies that there
exists a subsequence $(\boldsymbol{\tilde{v}_{k_{m}}})_{m}$ with
$\boldsymbol{\tilde{v}_{k_{m}}}\rightharpoonup\boldsymbol{v_{0,0}}\in V_{p}$
and this implies $\boldsymbol{\tilde{v}_{k_{m}}}\to\boldsymbol{v_{0,0}}$ in
$L^{p}(\Omega)^{N}$. This is a contradiction and yields the claim. ∎

######  Theorem 6.4 (Convergence for smooth solutions).

Let $\boldsymbol{v_{\mu_{0},\delta}}$ be the solution of
$G_{\mu_{0},\delta}(\boldsymbol{v})=0$ with the variables
$\mu_{0},\delta\in(0,\infty)$ or $(\mu_{0},\delta)=(0,0)$ as in the definition
of the operator $G$. Assume $\boldsymbol{v_{0,0}}\in V_{2}$. Then there exists
$\tilde{c}\in\mathbb{R}$ with

$\displaystyle\|\boldsymbol{v_{0,0}}-\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{2}}^{2}\leq\tilde{c}(|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}).$

###### Proof.

We follow the idea in [4, Theorem 4.1]. We set
$\boldsymbol{\tilde{v}}:=\boldsymbol{v_{0,0}}-\boldsymbol{v_{\mu_{0},\delta}}$
to shorten the notation. We use in the upcoming calculation the monotonicity
of $S^{p}$ and $S^{s}$,
$J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{\mu_{0},\delta}})=0$, and the main
theorem of calculus:

$\displaystyle\quad\frac{1}{2}\|\boldsymbol{\tilde{v}}\|_{V_{2}}^{2}$
$\displaystyle=\int_{0}^{1}\int_{\Omega}\frac{1}{r}|r\nabla\boldsymbol{\tilde{v}}|^{2}\,dx\,dr$
$\displaystyle\leq\int_{0}^{1}\int_{\Omega}\frac{1}{r}\big{(}S^{p}(D\boldsymbol{v_{\mu_{0},\delta}}+rD\boldsymbol{\tilde{v}})-S^{p}(D\boldsymbol{v_{\mu_{0},\delta}})\big{)}:(D\boldsymbol{v_{\mu_{0},\delta}}+rD\boldsymbol{\tilde{v}}-D\boldsymbol{v_{\mu_{0},\delta}})\,dx\,dr$
$\displaystyle\quad+\int_{0}^{1}\int_{\Gamma_{b}}\frac{1}{r}\big{(}S^{s}(\boldsymbol{v_{\mu_{0},\delta}}+r\boldsymbol{\tilde{v}})-S^{s}(\boldsymbol{v_{\mu_{0},\delta}})\big{)}\cdot(\boldsymbol{v_{\mu_{0},\delta}}+r\boldsymbol{\tilde{v}}-\boldsymbol{v_{\mu_{0},\delta}})\,ds\,dr$
$\displaystyle\quad+\int_{0}^{1}\int_{\Omega}\frac{1}{r}|r\nabla\boldsymbol{\tilde{v}}|^{2}\,dx\,dr$
$\displaystyle=\int_{0}^{1}\frac{1}{r}\big{(}J^{\prime}_{\mu_{0},\delta}(\boldsymbol{v_{\mu_{0},\delta}}+r\boldsymbol{\tilde{v}})-J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v_{\mu_{0},\delta}})\big{)}(\boldsymbol{v_{\mu_{0},\delta}}+r\boldsymbol{\tilde{v}}-\boldsymbol{v_{\mu_{0},\delta}})\,dr$
$\displaystyle=\int_{0}^{1}J^{\prime}_{\mu_{0},\delta}(\boldsymbol{v_{\mu_{0},\delta}}+r\boldsymbol{\tilde{v}})\boldsymbol{\tilde{v}}\,dr=J_{\mu_{0},\delta}(\boldsymbol{v_{0,0}})-J_{\mu_{0},\delta}(\boldsymbol{v_{\mu_{0},\delta}}).$

In the last step, we inserted $\boldsymbol{\tilde{v}}$. We use

$\displaystyle(|D\boldsymbol{v}|^{2}+\delta^{2})^{r/2}=\|(D\boldsymbol{v},\delta)\|_{2}^{r}\leq\|(D\boldsymbol{v},\delta)\|_{r}^{r}=|D\boldsymbol{v}|^{r}+\delta^{r}\quad\text{
for }r\in\\{s,p\\}$

to bound the functional $J_{\mu_{0},\delta}$ by $J_{\mu_{0},0}$ for all
$\boldsymbol{v}\in V_{2}$ with

$\displaystyle J_{\mu_{0},\delta}(\boldsymbol{v})$
$\displaystyle=\int_{\Omega}\frac{1}{p}(|D\boldsymbol{v}|^{2}+\delta^{2})^{p/2}\,dx+\int_{\Gamma_{b}}\frac{1}{s}(|\boldsymbol{v}|^{2}+\delta^{2})^{s/2}\,ds+\mu_{0}\|\boldsymbol{v}\|_{V_{2}}^{2}$
$\displaystyle\leq\frac{1}{p}\int_{\Omega}|D\boldsymbol{v}|^{p}+\delta^{p}\,dx+\frac{1}{s}\int_{\Gamma_{b}}|\boldsymbol{v}|^{s}+\delta^{s}\,ds+\mu_{0}\|\boldsymbol{v}\|_{V_{2}}^{2}$
$\displaystyle=J_{0,0}(\boldsymbol{v})+|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}\|\boldsymbol{v}\|_{V_{2}}^{2}.$

Hence, we obtain with the minimizer $\boldsymbol{v_{0,0}}$ of $J_{0,0}$ and
the definitions of $J_{0,0}$ and $J_{\mu_{0},\delta}$

$\displaystyle\quad
J_{\mu_{0},\delta}(\boldsymbol{v_{0,0}})-J_{\mu_{0},\delta}(\boldsymbol{v_{\mu_{0},\delta}})$
$\displaystyle\leq
J_{0,0}(\boldsymbol{v_{0,0}})+|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}\|\boldsymbol{v_{0,0}}\|_{V_{2}}^{2}-J_{\mu_{0},\delta}(\boldsymbol{v_{\mu_{0},\delta}})$
$\displaystyle\leq
J_{0,0}(\boldsymbol{v_{\mu_{0},\delta}})+|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}\|\boldsymbol{v_{0,0}}\|_{V_{2}}^{2}-J_{\mu_{0},\delta}(\boldsymbol{v_{\mu_{0},\delta}})$
$\displaystyle=\frac{1}{p}\int_{\Omega}|D\boldsymbol{v_{\mu_{0},\delta}}|^{p}-(\delta^{2}+|D\boldsymbol{v_{\mu_{0},\delta}}|^{2})^{p/2}\,dx+\frac{1}{s}\int_{\Gamma_{b}}|\boldsymbol{v_{\mu_{0},\delta}}|^{s}-(\delta^{2}+|\boldsymbol{v_{\mu_{0},\delta}}|^{2})^{s/2}\,ds$
$\displaystyle\quad+|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}\big{(}\|\boldsymbol{v_{0,0}}\|_{V_{2}}^{2}-\|\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{2}}^{2}\big{)}$
$\displaystyle\leq|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}\|\boldsymbol{v_{0,0}}\|_{V_{2}}^{2}.$

Thereby, we obtain with $c_{2}\in\mathbb{R}$

$\displaystyle\|\boldsymbol{v_{0,0}}-\boldsymbol{v_{\mu_{0},\delta}}\|_{V_{2}}^{2}\leq
c_{2}(|\Omega|\delta^{p}+|\Gamma_{b}|\delta^{s}+\mu_{0}).$

∎

##  7\. Numerical experiment

In this section, we describe the numerical experiment ISMIP-HOM $B$ formulated
in [25] and test it with our algorithms.

###  7.1. Computational domain

A unit square is too simple to represent real-world glaciers. Instead, in the
ISMIP-HOM experiments, see [25], a sinusoidal bedrock, see Fig. 1, is
suggested, and Dirichlet boundary conditions are imposed at the bedrock. We
set $\alpha:=0.5$ ∘, $L:=5000$ m, and define the upper and lower boundary of
the domain by $z_{s}:[0,L]\to\mathbb{R}$, $z_{b}:[0,L]\to\mathbb{R}$,

$\displaystyle z_{s}(x)=-\tan(\alpha)x\text{, }\quad\text{and
}z_{b}(x)=z_{s}(x)-1000+500\sin(\omega x)$

with $\omega=2\pi/L$. The experiment needs periodic boundary conditions on the
left and the right side of the domain and $\sigma\cdot\boldsymbol{n}=0$ at the
surface.

###  7.2. Parameters

In this subsection, we discuss setting constants and forcing boundary
conditions. We set the constants for the experiment corresponding to the
ISMIP-HOM experiments, see [25]: The ice parameter
$B:=0.5\cdot(10^{-16})^{-1/3}$ Pa-3a-1, the density $\rho:=910$ kg m-3, the
gravitational acceleration $g:=9.81$ m s-2, and the seconds per year
$31\,556\,926$. The nonlinear viscosity is
$(0.5|D\boldsymbol{v}|^{2}+\delta^{2})^{(p-2)/2}$ with $p:=4/3$. Periodic
boundary conditions have some difficulties in implementation for unstructured
grids. Moreover, they are not necessary for real-world applications. Thus, we
extend our domain by three copies to the left and the right. We apply
Dirichlet boundary conditions at these boundaries. This approach is suggested
in the Supplement of [25].

Figure 1. Considered domain in our experiment.

The resulting domain is shown in Fig. 1.

We rotate the gravity to obtain a horizontal upper surface. Formally, we
should rotate the left and right boundaries, too. However, this only changes
the position in $x$-direction by a maximum of $1500\textup{
m}\cdot\sin(0.5\textup{ \textdegree})\approx 13\textup{ m}$. This is
neglectable compared to the horizontal extent of $35\,000\textup{ m}$. Thus,
we simplify the domain with vertical left and right boundaries and a
horizontal surface boundary. Moreover, this approach is more flexible because
we can simply change the angle $\alpha$ without changing the domain to
generate different experiments. Furthermore, changes between using periodic
boundary conditions and the variant with copies of the glacier are easier. We
have eleven grid points in the vertical direction and $21$ grid points in the
horizontal direction for each copy of the glacier. Regarding applications, the
horizontal resolution must be smaller than the height of the glaciers. Our
resolution fulfills this condition because $5000/20=250<1000$.

The model is implemented in FEniCS, [26], version $2019.1.0$. FEniCS is a tool
that solves partial differential equations with finite elements. We use the
Taylor-Hood element $P_{2}-P_{1}$ for the velocity $\boldsymbol{v}$ and the
pressure $\pi$.

The value $\delta^{2}$ should be small compared to typical values of
$0.5|D\boldsymbol{v}|^{2}$. The range of $|D\boldsymbol{v}|$ is typically
between $10^{-7}$ /s and $10^{-11}$ /s. The calculations are done in years.
Hence, we obtain as a necessary condition

$\displaystyle\delta^{2}\leq 0.5(3.16\cdot
10^{-4}\textup{/a})^{2}\cdot\mathrm{eps}$

with the machine precision $\mathrm{eps}:=10^{-16}$. By rounding down to a
factor of $10$, we obtain $\delta:=10^{-12}\textup{/a}$. We should choose
$\mu_{0}$ such that

$\displaystyle\mu_{0}\leq(0.5|D\boldsymbol{v}|^{2}+\delta^{2})^{(p-2)/2}\cdot\mathrm{eps}$

is fulfilled. We obtain an upper bound by inserting the maximum typical value
of $|D\boldsymbol{v}|$ for $|D\boldsymbol{v}|$ and $\delta$. We conclude

$\displaystyle\mu_{0}\leq(2\cdot(10^{-7}\cdot
31\,556\,926\textup{/a})^{2})^{(4/3-2)/2}\cdot 10^{-16}\approx 3.69\cdot
10^{-17}\textup{a}^{2/3}.$

Consequently, we choose $\mu_{0}:=10^{-17}$ $\textup{kg
a/}(\textup{m}\textup{s}^{2})$. The unit for $\mu_{0}$ follows from
$[\mu_{0}|\nabla\boldsymbol{v}|^{2}]=[\rho\boldsymbol{g}\cdot\boldsymbol{v}]$.
To focus on the different solvers, we do not discuss different values of
$\delta$, even though [4] discussed faster convergence of Newton’s method and
higher accuracy dependent on the resolution for greater $\delta$ values.

###  7.3. Numerical solvers and computational effort

In this subsection, we describe the numerical algorithms to solve the problem
and the numerical effort for each iteration. The Picard iteration was
suggested for modeling glaciers in [27]. The Picard iteration is: For given
$\boldsymbol{v_{k}}\in V_{2}$ find $\boldsymbol{v_{k+1}}\in V_{2}$ with

$\displaystyle(B(|D\boldsymbol{v_{k}}|^{2}+\delta^{2})^{(p-2)/2}D\boldsymbol{v_{k+1}},\nabla\boldsymbol{\phi})$
$\displaystyle+(\tau(|\boldsymbol{v_{k}}|^{2}+\delta^{2})^{(s-2)/2}\boldsymbol{v_{k+1}},\boldsymbol{\phi})=-(\rho\boldsymbol{g},\boldsymbol{\phi})$

for all $\boldsymbol{\phi}\in V_{2}$. This algorithm is said to be often
globally convergent in practice but slow, [10]. It is, for example, used for
the Navier-Stokes equations with proved convergence theory, [16]. We will
compare the following algorithms:

  * •

The Picard iteration.

  * •

The Picard iteration with approximations of exact step sizes.

  * •

Newton’s method with approximations of exact step sizes.

  * •

Newton’s method with Armijo step sizes, see equation $(\ref{Armijo})$.

We can calculate the exact step size arbitrarily precise because we have a
convex functional. Let $0\leq a<b<\infty$. For the actual velocity
$\boldsymbol{v_{k}}$ and the direction $\boldsymbol{w_{k}}$, we can calculate

$\displaystyle\min_{t\in[a,b]}J_{help}(t):=\min_{t\in[a,b]}J_{\mu_{0},\delta}(\boldsymbol{v_{k}}+t\boldsymbol{w_{k}}).$

We set $b:=(a+b)/2$ for $J_{help}^{\prime}((a+b)/2)\geq 0$ and $a:=(a+b)/2$
for $J_{help}^{\prime}((a+b)/2)<0$, see [28, Algorithm 6.1], until we obtain
the wanted accuracy. Finally, we set $t:=(a+b)/2$. The chain rule yields for
$t\in(a,b)$

(20) $\displaystyle J_{help}^{\prime}(t)=\langle
G(\boldsymbol{v_{k}}+t\boldsymbol{w_{k}}),\boldsymbol{w_{k}}\rangle_{V_{2}^{*},V_{2}}.$

Furthermore, we can do these calculations without divergence-free elements:

######  Remark 7.1.

Only the calculation of the direction changes without divergence-free
elements. The step-size control is identical.

###### Proof.

We have to adapt the calculation of the Newton step and the Picard iteration,
respectively. After we do this, we have a divergence-free velocity
$\boldsymbol{v}$ and a divergence-free direction $\boldsymbol{w}$. Hence, we
have for the pressure $\pi$ and the direction $\psi$

$\displaystyle
J_{\mu_{0},\delta}^{\prime}(\boldsymbol{v})\boldsymbol{w}=\langle
G(\boldsymbol{v}),\boldsymbol{w}\rangle_{V_{2}^{*},V_{2}}=\langle
G(\boldsymbol{v}),\boldsymbol{w}\rangle_{V_{2}^{*},V_{2}}-(\pi,\mathrm{div}\boldsymbol{w})-(\mathrm{div}\boldsymbol{v},\psi).$

Thus, the claim follows, and minimizing $J_{\mu_{0},\delta}$ is sufficient. ∎

Now, we consider the computational effort. The computation consists of three
parts: Assembling the matrices, solving the linear system of equations, and
calculating the step size. The step size calculation consists of calculating
the value of an integral. Each integral calculation has an effort of $O(n)$
for the number of grid points $n$. In the implementation, the Armijo step size
needs a maximum of $20$ integral evaluations for each iteration. The
approximation of an exact step size needs $25$ integral evaluations for each
iteration. Solving the sparse stiffness matrix with the $LU$ decomposition has
an effort of $O(nm^{2})$ with the bandwidth $m$. We have for the bandwidth as
a lower bound the number of grid points in vertical direction. Thus, we have
$m\geq 11$, and the step size calculation is negligible compared to solving
the problem.

###  7.4. Results and interpretation

The ISMIP-HOM $B$ experiment measures the norms of the surface velocity
$v_{r}$:

$\displaystyle
v_{r}:=\sqrt{(v_{1}\cdot\cos(\alpha))^{2}-(v_{2}\cdot\sin(\alpha))^{2}}.$

The surface velocities are approximately the same for all our methods, see
Fig. 2. Moreover, the velocities are similar to [25, Figure 6].

Figure 2. The velocity norms at the surface are shown in the left figure. The
velocity field is shown in the right figure.

Now, we concentrate on our main purpose, increasing the convergence speed. All
step size controls used the convex functional, see Lemma 5.2, as the
minimization term.

Figure 3. The relative residual norms are visualized in the left figure, and
the functional values are visualized in the right figure.

All methods reduce the residual error compared to the initial residual error
by approximately $10^{3}$ or more. The residual error of the Picard iteration
decreases linearly.

Interestingly, Newton’s method with Armijo step sizes, Newton’s method with
approximations of exact step sizes, and the Picard iteration with
approximations of exact step sizes behave similarly for the first few
iterations. However, the reduction of the norm stops earlier for Newton’s
method with Armijo step sizes than for all other methods. It seems like the
functional for step size control is unsuitable for reducing the error after
some iterations because then the functional values are nearly constant.

We explain this behavior: An unique solution $\boldsymbol{v}$ exists for the
regularized $p$-Stokes equations. One unique solution $\boldsymbol{v_{h}}$
should exist for the discretization $G_{h}$ of $G$ with
$G_{h}(\boldsymbol{v_{h}})=0$. Also, a unique minimum
$\boldsymbol{\tilde{v}_{h}}$ should exist for the discretized convex
functional $J_{h}$ of $J$. Because of $J_{h}\to J$ and $G_{h}\to G$ for the
mesh size $h\to 0$, we have

$\displaystyle\boldsymbol{v_{h}}\to\boldsymbol{v}\leftarrow\boldsymbol{\tilde{v}_{h}}\quad\text{
for }h\to 0.$

But we can have $\boldsymbol{v_{h}}\neq\boldsymbol{\tilde{v}_{h}}$ for $h>0$.

######  Remark 7.2.

Let $f:\mathbb{R}\to\mathbb{R}$, $f(x)=x^{2}$. We interpolate this function
with linear elements as well as the derivative. Let the interpolation points
be $-0.5$ and $1$. That yields for the discretized function $f_{h}$:

$\displaystyle f_{h}(-0.5)=\min_{x\in\mathbb{R}}f_{h}(x)\text{, }\quad
f^{\prime}_{h}(0)=\min_{x\in\mathbb{R}}|f^{\prime}_{h}(x)|.$

Thus, the solution that we want to find is at different points for the
function and the derivative, respectively.

Hence, the quality of Newton’s method should be dependent on the mesh
resolution if we estimate the quality with the residual norm. We consider the
residual norm because it is of more relevance in practice as minimizing the
functional.

Figure 4. Relative residual norms for different mesh resolutions. The length
divided by the number is the resolution in this direction. The first number
describes the horizontal direction, and the second number the vertical
direction.

Newton’s method with Armijo step sizes depends on the mesh resolution, see
Fig. 4. Thus, our assumption is justified. The approximation of exact step
sizes seems less dependent on the mesh resolution. That is also suggested by
the formulation of the directional derivative of the convex functional in
equation $(\ref{DerivativeHelpFunc})$.

##  8\. Conclusion

In conclusion, it is possible to obtain a global convergent Newton method with
step size control with a convex functional. Moreover, the exact step size
control is also good. More experiments should be done in three dimensions,
with nonlinear sliding and the interaction of the glacier with water.
Furthermore, an implementation in ice sheet models could be interesting.

### Code availability

The source code for the simulation can be found at
https://doi.org/10.5281/zenodo.7805479.

## Acknowledgments

I thank my supervisor Prof. Thomas Slawig from Kiel University, for many
helpful discussions about mathematical and glaciological problems. Moreover, I
thank Prof. Angelika Humbert and Dr. Thomas Kleiner from Alfred-Wegener-
Institut in Bremerhaven, Dr. Martin Rückamp from the Bavarian Academy of
Sciences and Humanities, and Prof. Andreas Rademacher from the University of
Bremen for helpful discussions.

## References

  * [1] Fischler, Y., Rückamp, M., Bischof, C., Aizinger, V., Morlighem, M., and Humbert, A.: A scalability study of the Ice-sheet and Sea-level System Model (ISSM, version 4.18).  _Geoscientific Model Development_ 15, 3753–3771 (2022)
  * [2] Chen, Q., Gunzburger, M., and Perego, M.: Well-Posedness Results for a Nonlinear Stokes Problem Arising in Glaciology.  _SIAM Journal on Mathematical Analysis_ 45, 2710–2733 (2013)
  * [3] Jouvet, G. and Rappaz, J.: Analysis and Finite Element Approximation of a Nonlinear Stationary Stokes Problem Arising in Glaciology.  _Advances in Numerical Analysis_ 2011, 1–24 (2011)
  * [4] Hirn, A.: Finite element approximation of singular power-law systems.  _Mathematics of Computation_ 82, 1247–1268 (2013)
  * [5] de Diego, G. G., Farrell, P. E., and Hewitt, I. J.: On the Finite Element Approximation of a Semicoercive Stokes Variational Inequality Arising in Glaciology.  _SIAM Journal on Numerical Analysis_ 61, 1–25 (2023)
  * [6] Fowler, A. and Ng, F., eds.: _Glaciers and Ice Sheets in the Climate System_.  Springer International Publishing (2021)
  * [7] Casas, E. and Fernandez, L.: Distributed control of systems governed by a general class of quasilinear elliptic equations.  _J of Diff Equations_ 104, 20–47 (1993)
  * [8] Arada, N.: Distributed control for multistate modified Navier-Stokes equations.  _ESAIM: Control, Optimisation and Calculus of Variations_ 19, 219–238 (2012)
  * [9] Arada, N.: Optimal Control of Shear-Thinning Fluids.  _SIAM Journal on Control and Optimization_ 50, 2515–2542 (2012)
  * [10] Fraters, M. R. T., Bangerth, W., Thieulot, C., Glerum, A. C., and Spakman, W.: Efficient and practical Newton solvers for non-linear Stokes systems in geodynamic problems.  _Geophysical Journal International_ 218, 873–894 (2019)
  * [11] Allen, J., Leibs, C., Manteuffel, T., and Rajaram, H.: A Fluidity-Based First-Order System Least-Squares Method for Ice Sheets.  _SIAM Journal on Scientific Computing_ 39, B352–B374 (2017)
  * [12] Bons, P. D., Kleiner, T., Llorens, M.-G., Prior, D. J., Sachau, T., Weikusat, I., and Jansen, D.: Greenland Ice Sheet: Higher Nonlinearity of Ice Flow Significantly Reduces Estimated Basal Motion.  _Geophysical Research Letters_ 45, 6542–6548 (2018)
  * [13] MacGregor, J. A., Fahnestock, M. A., Catania, G. A., Aschwanden, A., Clow, G. D., Colgan, W. T., Gogineni, S. P., Morlighem, M., Nowicki, S. M. J., Paden, J. D., Price, S. F., and Seroussi, H.: A synthesis of the basal thermal state of the Greenland Ice Sheet.  _Journal of Geophysical Research: Earth Surface_ 121, 1328–1350 (2016)
  * [14] Chipot, M.: _$\ell$ Goes to Plus Infinity_.  Birkhäuser Basel (2002)
  * [15] Amrouche, C. S. and Girault, V.: Decomposition of vector spaces and application to the Stokes problem in arbitrary dimension.  _Czechoslovak Mathematical Journal_ 44, 109–140 (1994)
  * [16] Girault, V. and Raviart, P.-A.: _Finite Element Methods for Navier-Stokes Equations_.  Springer Series in Computational Mathematics 5, New York (1986)
  * [17] Browder, F. E.: Nonlinear elliptic boundary value problems.  _Bulletin of the American Mathematical Society_ 69, 862 – 874 (1963)
  * [18] Diening, L., Ebmeyer, C., and Růžička, M.: Optimal Convergence for the Implicit Space-Time Discretization of Parabolic Systems with p-Structure.  _SIAM Journal on Numerical Analysis_ 45, 457–472 (2007)
  * [19] Hinze, M., Pinnau, R., Ulbrich, M., and Ulbrich, S.: _Optimization with PDE Constraints_.  Springer Netherlands (2009)
  * [20] Málek, J., Nečas, J., Rokyta, M., and Ružička, M.: _Weak and Measure-valued Solutions to Evolutionary PDEs_.  Chapman & Hall (1996)
  * [21] Polyak, R. A.: Regularized Newton method for unconstrained convex optimization.  _Mathematical Programming_ 120, 125–145 (2007)
  * [22] Bertsekas, D. M.: _Constrained Optimization and Lagrange Multiplier Methods_.  Academic Press, Paris (1982)
  * [23] Haase, M.: _Functional analysis an elementary introduction_.  American Mathematical Society (2014).  ISBN 9780821891711
  * [24] Zeidler, E.: _Nonlinear functional analysis and its application II/B: Nonlinear monotone operators_.  Springer, New York (1990)
  * [25] Pattyn, F., Perichon, L., Aschwanden, A., Breuer, B., de Smedt, B., Gagliardini, O., Gudmundsson, G. H., Hindmarsh, R. C. A., Hubbard, A., Johnson, J. V., Kleiner, T., Konovalov, Y., Martin, C., Payne, A. J., Pollard, D., Price, S., Rückamp, M., Saito, F., Souček, O., Sugiyama, S., and Zwinger, T.: Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP–HOM).  _The Cryosphere_ 2, 95–108 (2008)
  * [26] Logg, A., Mardal, K.-A., and Wells, G., eds.: _Automated Solution of Differential Equations by the Finite Element Method_.  Springer Berlin Heidelberg (2012)
  * [27] Colinge, J. and Rappaz, J.: A strongly nonlinear problem arising in glaciology.  _ESAIM: Mathematical Modelling and Numerical Analysis_ 33, 395–406 (1999)
  * [28] Schölkopf, B. and Smola, A. J.: _Learning with Kernels : Support Vector Machines, Regularization, Optimization, and Beyond_.  MIT Press (2001)