# Guaranteed, locally efficient, and robust a posteriori estimates for nonlinear elliptic problems in iteration-dependent norms. An orthogonal decomposition result based on iterative linearization

Koondanibha Mitra, Martin Vohralík

► To cite this version:

Koondanibha Mitra, Martin Vohralík. Guaranteed, locally efficient, and robust a posteriori estimates for nonlinear elliptic problems in iteration-dependent norms. An orthogonal decomposition result based on iterative linearization. 2023. hal-04156711

HAL Id: hal-04156711

https://inria.hal.science/hal-04156711

Preprint submitted on 9 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Guaranteed, locally efficient, and robust a posteriori estimates for nonlinear elliptic problems in iteration-dependent norms. An orthogonal decomposition result based on iterative linearization

K. Mitra¹ and M. Vohralík²,³

¹Hasselt University, Belgium

²Inria, 2 rue Simone Iff, 75589 Paris, France

³CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France

July 9, 2023

Abstract

We consider numerical approximations of nonlinear, monotone, and Lipschitz-continuous elliptic problems, with gradient-dependent or gradient-independent diffusivity. For this purpose, we first employ a symmetric iterative linearization, with examples including respectively the Kačanov, Newton, Zarantonello and fixed-point (Picard), *L*, and *M*-schemes. We design an iteration-dependent norm, stemming from the concrete linearization and finite element discretization, which leads to an orthogonal decomposition of the total error into the linearization and discretization components. This norm is of reaction-diffusion type, so that available results from a posteriori analysis of linear (singularly perturbed) reaction-diffusion problems enable us to obtain guaranteed and locally efficient error bounds robust with respect to the strength of the nonlinearity. Moreover, the clear distinction and orthogonality of the error components leads to a design of adaptive stopping criteria for iterative linearizations ensuring a natural balance between discretization and linearization components. Under appropriate assumptions, we also show that the iteration-dependent norm is, up to a multiplicative factor tending to one, equivalent to the iteration-independent natural discrete energy norm. In an appendix, we address the Newton (nonsymmetric) linearization for the gradient-independent diffusivity case. Numerical experiments illustrate the theoretical developments.

Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>1.1</td><td>H<sub>0</sub><sup>1</sup>(&Omega;)-error-norm and energy difference a posteriori estimates</td><td>3</td></tr><tr><td>1.2</td><td>Dual norm of the residual a posteriori estimates</td><td>4</td></tr><tr><td>1.3</td><td>Our goals</td><td>4</td></tr><tr><td>1.4</td><td>Motivation</td><td>4</td></tr><tr><td>1.4.1</td><td>Energy norm estimates in the linear case</td><td>4</td></tr><tr><td>1.4.2</td><td>A first attempt in the nonlinear case</td><td>5</td></tr><tr><td>1.5</td><td>The pathway on an example</td><td>5</td></tr><tr><td>1.5.1</td><td>Iterative linearization</td><td>5</td></tr><tr><td>1.5.2</td><td>Iteration-dependent energy norm</td><td>6</td></tr><tr><td>1.5.3</td><td>Residuals</td><td>6</td></tr><tr><td>1.5.4</td><td>Orthogonal error decomposition</td><td>6</td></tr><tr><td>1.5.5</td><td>A posteriori error estimates</td><td>6</td></tr><tr><td>1.5.6</td><td>Balancing of error components and stopping criteria</td><td>7</td></tr><tr><td>1.6</td><td>Structure of the present contribution</td><td>7</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>2</td>
      <td>Basic mathematical tools and notation</td>
      <td>7</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Two classes of nonlinear elliptic problems</td>
      <td>8</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Gradient-dependent diffusivity</td>
      <td>8</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Gradient-independent diffusivity</td>
      <td>9</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Finite element discretization and iterative linearization</td>
      <td>10</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Finite element discretization</td>
      <td>10</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Iteration-dependent scalar product</td>
      <td>10</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Iterative linearization</td>
      <td>11</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Continuous-level counterpart of (4.7)</td>
      <td>12</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Iteration-dependent norm</td>
      <td>12</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Examples of common linearization schemes</td>
      <td>12</td>
    </tr>
    <tr>
      <td>4.6.1</td>
      <td>Gradient-dependent diffusivity</td>
      <td>13</td>
    </tr>
    <tr>
      <td>4.6.2</td>
      <td>Gradient-independent diffusivity</td>
      <td>13</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Main results</td>
      <td>14</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>An orthogonal decomposition of the error</td>
      <td>14</td>
    </tr>
    <tr>
      <td>5.1.1</td>
      <td>Discretization residual and discretization error</td>
      <td>14</td>
    </tr>
    <tr>
      <td>5.1.2</td>
      <td>Total residual and total error</td>
      <td>15</td>
    </tr>
    <tr>
      <td>5.1.3</td>
      <td>An orthogonal decomposition of the total residual/error</td>
      <td>15</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Guaranteed and locally efficient a posteriori error bounds robust with respect to the strength of the nonlinearity</td>
      <td>16</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Balancing of the linearization and discretization components and stopping criteria for iterative linearizations</td>
      <td>19</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Numerical experiments</td>
      <td>19</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Gradient-dependent diffusivity</td>
      <td>20</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Gradient-independent diffusivity: the Richards equation</td>
      <td>22</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Well-posedness of the models (3.1) and (3.3)</td>
      <td>27</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Well-posedness of the gradient-dependent diffusivity model (3.1)</td>
      <td>27</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Well-posedness of the gradient-independent diffusivity model (3.3)</td>
      <td>28</td>
    </tr>
    <tr>
      <td>8</td>
      <td>A posteriori error estimates of Theorem 5.7</td>
      <td>29</td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>Equilibrated flux and potential of Definition 5.4, proof of Lemma 5.6</td>
      <td>30</td>
    </tr>
    <tr>
      <td>8.1.1</td>
      <td>Orthogonal projections</td>
      <td>30</td>
    </tr>
    <tr>
      <td>8.1.2</td>
      <td>Trace and inverse inequalities</td>
      <td>30</td>
    </tr>
    <tr>
      <td>8.1.3</td>
      <td>Weights</td>
      <td>30</td>
    </tr>
    <tr>
      <td>8.1.4</td>
      <td>Shorhand notations</td>
      <td>31</td>
    </tr>
    <tr>
      <td>8.1.5</td>
      <td>Definition 5.4</td>
      <td>31</td>
    </tr>
    <tr>
      <td>8.1.6</td>
      <td>Proof of Lemma 5.6</td>
      <td>31</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Proof of global reliability estimate in Theorem 5.7</td>
      <td>32</td>
    </tr>
    <tr>
      <td>8.3</td>
      <td>Local efficiency bounds</td>
      <td>33</td>
    </tr>
    <tr>
      <td>8.4</td>
      <td>Proof of efficiency estimates in Theorem 5.7</td>
      <td>36</td>
    </tr>
    <tr>
      <td colspan="2">Appendices</td>
      <td>37</td>
    </tr>
    <tr>
      <td colspan="2">A Iterative linearizations in abstract spaces: well-posedness and consistency</td>
      <td>37</td>
    </tr>
    <tr>
      <td colspan="2">B Error estimates in a fixed norm</td>
      <td>39</td>
    </tr>
  </tbody>
</table>
## 1.2 Dual norm of the residual a posteriori estimates

To circumvent the dependence on $\lambda_M/\lambda_m$, observe that

$$ \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(u_\ell) - \mathcal{R}(u), \varphi \rangle}{\|\nabla \varphi\|} = \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(u_\ell), \varphi \rangle}{\|\nabla \varphi\|} = \|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)}, \quad (1.6) $$

since $u \in H_0^1(\Omega)$ solves (1.1). One can then consider $\|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)}$ directly as an error measure, since it in particular vanishes if and only if $u_\ell = u$ (it is equivalent to $\text{dist}(u_\ell, u)$ by (1.2)). Moreover, $\|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)}$, though defined as a global dual norm, is localizable by the orthogonality to the lowest-order shape functions from (1.3); this is explicitly recalled in, e.g., [6, 11, 16], see also the references therein. Furthermore, it can be a posteriori estimated robustly by the computable estimator $\eta(u_\ell)$, i.e.,

$$ \|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)} \le \eta(u_\ell) \le C \|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)}, \quad (1.7) $$

where the constant is as in (1.4) independent of the global Lipschitz/monotonicity ratio $\lambda_M/\lambda_m$, see [14, 15, 23, 24]. Nevertheless, $\|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)}$ might be too weak an error measure since it does not carry any physical units/scaling/energy information and may not correspond to an error metric either locally or globally.

## 1.3 Our goals

In the present contribution, we are interested in deriving a posteriori error estimates for finite element approximations (1.3) of (1.1) independent of the factor $\lambda_M/\lambda_m$, i.e., robust with respect to the strength of the nonlinearity, in contrast to (1.4) or (1.5). In improvement of (1.7), energy information and local norm correspondence are encoded. Based on iterative linearization, we also derive an orthogonal decomposition result of the total error into the discretization and linearization error components. This finally allows us to design an adaptive stopping criterion for iterative linearizations ensuring a natural balance between these error components.

## 1.4 Motivation

We start by motivating our approach.

### 1.4.1 Energy norm estimates in the linear case

To gain insight, let us turn our attention for a while to a linear problem. Consider the linear diffusion equation

$$ \langle \mathcal{R}(u), \varphi \rangle := (\bar{\mathcal{K}}\nabla u, \nabla \varphi) + \langle f, \varphi \rangle = 0 \quad \forall \varphi \in H_0^1(\Omega), \quad (1.8) $$

where $f \in H^{-1}(\Omega)$ and, for a.e. $x \in \Omega$, $\bar{\mathcal{K}}(x) \in \mathbb{R}^{d \times d}$ is symmetric and satisfies the uniform ellipticity and boundedness condition, i.e., there exist real constants $\lambda_M \ge \lambda_m > 0$ such that $\lambda_m |y|^2 \le y^T \bar{\mathcal{K}}(x)y \le \lambda_M |y|^2$ for all $y \in \mathbb{R}^d$. In this case, by Céa's lemma [17], one has for the finite element approximation $u_\ell \in V_\ell$ given by (1.3) the a priori characterization

$$ \|\nabla(u - u_\ell)\| \le \frac{\lambda_M}{\lambda_m} \|\nabla(u - \varphi_\ell)\| \quad \forall \varphi_\ell \in V_\ell. \quad (1.9) $$

The estimate above implies that $u_\ell$ is only a quasi-optimal approximation of $u$ from (1.8) in $V_\ell$ with respect to the intrinsic $H_0^1(\Omega)$-norm, not strictly leading to error minimization and involving the ratio $\lambda_M/\lambda_m$. Congruently, similar to (1.2) and (1.4), one only obtains the a posteriori characterizations

$$ \lambda_m \|\nabla(u - u_\ell)\| \le \|\mathcal{R}(u_\ell)\|_{H^{-1}(\Omega)} \le \lambda_M \|\nabla(u - u_\ell)\|, \quad (1.10a) $$

$$ \lambda_m \|\nabla(u - u_\ell)\| \le \eta(u_\ell) \le C \lambda_M \|\nabla(u - u_\ell)\|, \quad (1.10b) $$
which involve the ratio $\lambda_M/\lambda_m$ and are in particular not robust with respect to the inhomogeneity and anisotropy in $\bar{K}$.

In contrast to the above, introducing the energy norm $\|\cdot\|_{1,\bar{K}} := (\bar{K}\nabla\varphi, \nabla\varphi)^{\frac{1}{2}}$, which is arguably the natural norm for this system, one has by Céa's lemma that (compare with (1.9))

$$ \|u - u_{\ell}\|_{1, \bar{K}} = \min_{\varphi_{\ell} \in V_{\ell}} \|u - \varphi_{\ell}\|_{1, \bar{K}}. \qquad (1.11) $$

This entails that $u_{\ell} \in V_{\ell}$ is the optimal approximation of $u$ from (1.8) in $V_{\ell}$ with respect to the norm $\|\cdot\|_{1,\bar{K}}$; this equality crucially does not involve the ratio $\lambda_M/\lambda_m$ and is consequently robust with respect to $\lambda_M/\lambda_m$. Moreover, (1.11) motivates to introduce the energy dual norm

$$ \|\zeta\|_{-1, \bar{K}} := \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \zeta, \varphi \rangle}{\|\varphi\|_{1, \bar{K}}} $$

for $\zeta \in H^{-1}(\Omega)$ and consider $\|\mathcal{R}(u_{\ell})\|_{-1, \bar{K}}$ as an error measure alternative to $\|\mathcal{R}(u_{\ell})\|_{H^{-1}(\Omega)}$. In fact,

$$ \begin{aligned} \|\mathcal{R}(u_{\ell})\|_{-1, \bar{K}} &= \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(u_{\ell}), \varphi \rangle}{\|\varphi\|_{1, \bar{K}}} = \sup_{\varphi \in H_0^1(\Omega)} \frac{(\bar{K}\nabla u_{\ell}, \nabla\varphi) + \langle f, \varphi \rangle}{\|\varphi\|_{1, \bar{K}}} \\ &= \sup_{\varphi \in H_0^1(\Omega)} \frac{(\bar{K}\nabla(u_{\ell} - u), \nabla\varphi)}{\|\varphi\|_{1, \bar{K}}} = \|u - u_{\ell}\|_{1, \bar{K}}, \end{aligned} $$

so that (1.10a) is also turned into equality. Just as $\|u - u_{\ell}\|_{1, \bar{K}}$ can be estimated robustly a priori using (1.11), robust a posteriori error estimates of the form (compare with (1.10b))

$$ \|u - u_{\ell}\|_{1, \bar{K}} \le \bar{\eta}(u_{\ell}) \le C \|u - u_{\ell}\|_{1, \bar{K}} $$

become available at least under some conditions on $\bar{K}$, see [5, 12, 65] and the references therein. The use of the energy norm and the dual energy norm of the residual also lead to robust a posteriori estimates for singularly perturbed reaction-diffusion problems, see [1, 2, 42, 43, 60, 64] and the references therein.

## 1.4.2 A first attempt in the nonlinear case

The above approach cannot immediately be generalized to the nonlinear case. Indeed, let us consider a nonlinear extension of (1.8)

$$ \langle \mathcal{R}(u), \varphi \rangle := (\bar{K} \mathcal{D}(u) \nabla u, \nabla \varphi) + \langle f, \varphi \rangle = 0 \quad \forall \varphi \in H_0^1(\Omega), \qquad (1.12) $$

where $\mathcal{D}: \mathbb{R} \to (0, \infty)$ is a non-constant function. In this case, the norms $\|\cdot\|_{\pm 1, \bar{K}\mathcal{D}(u)}$ can no longer be directly computed since the exact solution $u \in H_0^1(\Omega)$ is unknown, and a modification of the approach from the linear case is necessary.

# 1.5 The pathway on an example

We now stick to the example (1.12) and describe our approach in this model setting.

## 1.5.1 Iterative linearization

Iterative linearizations are commonly used to approximate the nonlinear problems (1.1) and (1.3). One of the most straightforward is the fixed point (Picard) linearization. For the specific case (1.12), the finite element discretization combined with the Picard linearization reads: let
$u_{\ell}^i \in V_\ell$ be an approximation of $u_\ell$ at some linearization iteration $i \ge 0$. We then look for $u_{\ell}^{i+1} \in V_\ell$ such that

$$ (\bar{K}D(u_{\ell}^{i})\nabla u_{\ell}^{i+1}, \nabla \varphi_{\ell}) + \langle f, \varphi_{\ell} \rangle = 0 \quad \forall \varphi_{\ell} \in V_{\ell}. \qquad (1.13) $$

We now observe that $u_{\ell}^{i+1}$ is the finite element approximation of the continuous-level solution $u_{\langle \ell \rangle}^{i+1} \in H_0^1(\Omega)$ of the linear problem

$$ (\bar{K}D(u_{\ell}^{i})\nabla u_{\ell}^{i+1}, \nabla \varphi) + \langle f, \varphi \rangle = 0 \quad \forall \varphi \in H_{0}^{1}(\Omega). \qquad (1.14) $$

### 1.5.2 Iteration-dependent energy norm

Since $u_{\ell}^i \in V_\ell$ is known, at each iteration step $i \ge 0$, similarly to the linear case of Section 1.4.1, one can introduce the iteration-dependent energy norm $\|\cdot\|_{1,u_{\ell}^i}$ and the corresponding dual norm $\|\cdot\|_{-1,u_{\ell}^i}$ by

$$ \begin{aligned} \|\varphi\|_{1, u_{\ell}^i} &:= (\bar{K}D(u_{\ell}^i)\nabla\varphi, \nabla\varphi)^{\frac{1}{2}}, & \varphi \in H_0^1(\Omega), \\ \|\zeta\|_{-1, u_{\ell}^i} &:= \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \zeta, \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^i}}, & \zeta \in H^{-1}(\Omega). \end{aligned} $$

### 1.5.3 Residuals

Let the discretization residual corresponding to the $(i+1)$ linearization iteration be given by

$$ \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1}), \varphi \rangle := (\bar{K}D(u_{\ell}^{i})\nabla u_{\ell}^{i+1}, \nabla \varphi) + \langle f, \varphi \rangle, \quad \varphi \in H_{0}^{1}(\Omega). $$

Note that from (1.14), $\langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle$ uniformly vanishes if and only if there is no discretization error in the sense that $u_{\ell}^{i+1} = u_{\langle \ell \rangle}^{i+1}$. Using (1.12) and (1.14), one observes that, for all $\varphi \in H_0^1(\Omega)$, the total residual writes

$$ \langle \mathcal{R}(u_{\ell}^{i}), \varphi \rangle = (\bar{K}D(u_{\ell}^{i})\nabla u_{\ell}^{i}, \nabla \varphi) + \langle f, \varphi \rangle = (\bar{K}D(u_{\ell}^{i})\nabla(u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}), \nabla \varphi) + \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1}), \varphi \rangle. $$

### 1.5.4 Orthogonal error decomposition

Now, using the Galerkin orthogonality from (1.13)-(1.14) reveals the orthogonal decomposition

$$ \underbrace{\|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2}}_{\text{total residual/error}} = \underbrace{\|u_{\ell}^{i} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}}^{2}}_{\text{linearization error}} + \underbrace{\|\mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})\|_{-1, u_{\ell}^{i}}^{2}}_{\substack{\text{discretization residual/} \\ \text{error}}} \qquad (1.15) $$

A detailed proof of this result is given in Section 5.1.3 for the general case. Clearly, the first term on the right-hand-side is the linearization error measured by the difference of two consecutive iterates, whereas the second one represents the discretization error on step $(i+1)$ of the linearization. Consequently, $\|\mathcal{R}(u_{\ell}^{i})\|_{-1,u_{\ell}^{i}} = \|\|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}\|\|_{1,u_{\ell}^{i}}$ seems as the ideal error metric for the nonlinear problem (1.12), with orthogonal decomposition into the two error components.

### 1.5.5 A posteriori error estimates

The linearization error in (1.15) is directly computable. In addition, crucially, $\|\mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})\|_{-1,u_{\ell}^{i}}$ which equals $\|u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1}\|_{1,u_{\ell}^{i}}$ and stems from a linear (here diffusion and below possibly singularly perturbed reaction-diffusion) problem may be robustly estimated using the energy setting of the linear case of Section 1.4.1. Moreover, the total error $\|\mathcal{R}(u_{\ell}^{i})\|_{-1,u_{\ell}^{i}}$ converges to
$\|\mathcal{R}(u_\ell)\|_{-1, \bar{K}_D(u_\ell)}$ from Section 1.4.2 if the linearization is converging. This is shown for the general case in Appendix B.

### 1.5.6 Balancing of error components and stopping criteria

Finally, the orthogonal decomposition (1.15) allows us to naturally balance the two error components and design an adaptive stopping criterion for the linearization iterations. The linearization error component arising from our analysis has also been used in [61] to switch between linearization schemes, L-scheme and the Newton scheme to be more precise, when solving the time-discretized Richards equation.

## 1.6 Structure of the present contribution

In this manuscript, we consider two rather broad classes of nonlinear elliptic problems of the form (1.1); details are elaborated in Section 3 after we fix the setting in Section 2. We then in Section 4 describe the class of linearization schemes that we consider, including respectively the Kačanov, Newton, Zarantonello and the fixed-point (Picard), $L$, and $M$-schemes. Importantly, they all give rise to a symmetric linearization just as in (1.13) (of reaction-diffusion type in general). Our main results are collected in Section 5. First, an orthogonal decomposition result of the form (1.15) is shown in Section 5.1. Then, in Section 5.2, guaranteed and locally efficient a posteriori error estimates not featuring the global Lipschitz/monotonicity ratio $\lambda_M/\lambda_m$ are stated. Finally, an adaptive stopping criterion for iterative linearizations ensuring a natural balance between discretization and linearization components is designed in Section 5.3. In Section 6, we then present numerical results for two test problems and five different linearization schemes. They clearly illustrate the robustness as well as effectiveness of the derived a posteriori estimates. Next, in Section 7 functional analytic aspects of Section 3 are discussed in detail, and, subsequently, the details and the proof of the a posteriori estimates are presented in Section 8. Appendix A summarizes the well-posedness and consistency of the iterative linearizations introduced in Section 4 on an abstract level. The transition of our results from iteration-dependent norms to fixed norms is then discussed in Appendix B. Finally, Appendix C extends our results to the Newton linearization for the gradient-independent diffusivity case, which, in contrast to the above, gives rise to a linear but nonsymmetric (advection-reaction-diffusion) problem on each linearization step. The framework has to be extended but, under appropriate (smallness) assumptions, similar structural results as in the main body of the manuscript are obtained, with convincing numerical illustrations.

# 2 Basic mathematical tools and notation

Let $L^q(\Omega)$ be the space of functions that are Lebesgue-integrable with power $1 \le q \le \infty$. For $v \in L^q(\Omega)$, we denote the corresponding norm by $\|v\|_{L^q(\Omega)}$ and reduce this notation to $\|v\|$ in the specific case $q = 2$; there, the corresponding scalar product is denoted by $(v, w)$. Similar notation is used for vector-valued functions. On a subdomain $\omega \subseteq \Omega$, we use the notation $\|v\|_{\omega}$ and $(v, w)_{\omega}$. We denote by $W^{1,q}(\Omega)$ the Sobolev space of $L^q(\Omega)$ functions having a weak derivative also in $L^q(\Omega)$, with $H^1(\Omega) := W^{1,2}(\Omega)$; $W_0^{1,q}(\Omega)$ and $H_0^1(\Omega)$ are then the trace-free subspaces.

For an open Lipschitz domain $\omega \subseteq \Omega$ with diameter $h_\omega > 0$, we will repeatedly use the Poincaré inequality: for functions $u \in H^1(\omega)$ vanishing in a trace sense in a fixed relatively open subset of $\partial\omega$ with non-zero measure, or for functions $u \in H^1(\omega)$ satisfying $\int_\omega u = 0$ (zero-mean-value property), there holds for a constant $C_\omega > 0$

$$ \|u\|_{\omega} \le C_{\omega} h_{\omega} \|\nabla u\|_{\omega}. \qquad (2.1) $$
If $\omega$ is convex, then $C_\omega \le \pi^{-1}$ in the zero-mean-value case. The discussion of the other cases can be found in, e.g., [6, Section 2.1] and the references therein. We will further use the Young inequality which states that for $\rho \in \mathbb{R}$, $\rho > 0$, and $a, b \in \mathbb{R}$,

$$ab \le \frac{1}{2\rho}a^2 + \frac{\rho}{2}b^2. \qquad (2.2)$$

# 3 Two classes of nonlinear elliptic problems

For our analysis, we divide elliptic problems represented by (1.1) into two different classes.

## 3.1 Gradient-dependent diffusivity

In this case, the action of the nonlinear operator $\mathcal{R}$ applied on a function $\xi \in H_0^1(\Omega)$ is prescribed by

$$\langle \mathcal{R}(\xi), \varphi \rangle := (\sigma(\cdot, \nabla\xi), \nabla\varphi) + \langle f(\cdot, \xi), \varphi \rangle, \quad \varphi \in H_0^1(\Omega). \qquad (3.1)$$

We suppose the following to ensure the monotonicity and Lipschitz continuity of the operator $\mathcal{R}$ in the sense of (1.2):

**Assumption 3.1 (Gradient-dependent diffusivity).**

(A1) The flux function $\sigma: \Omega \times \mathbb{R}^d \to \mathbb{R}^d$ is bounded with respect to its first argument $x \in \Omega$ for a fixed second argument $y \in \mathbb{R}^d$. Moreover, it satisfies for real constants $\sigma_M \ge \sigma_m > 0$ the Lipschitz-continuity condition

$$|\sigma(x, y) - \sigma(x, z)| \le \sigma_M |y - z| \quad \text{for a.e. } x \in \Omega \text{ and all } y, z \in \mathbb{R}^d$$

and the monotonicity condition

$$(\sigma(x, y) - \sigma(x, z)) \cdot (y - z) \ge \sigma_m |y - z|^2 \quad \text{for a.e. } x \in \Omega \text{ and all } y, z \in \mathbb{R}^d.$$

(A2) The reaction/source function $f: \Omega \times \mathbb{R} \to \mathbb{R}$ is bounded with respect to its first argument $x \in \Omega$ for a fixed second argument $\xi \in \mathbb{R}$. It is Lipschitz continuous and increasing with respect to its second argument $\xi \in \mathbb{R}$, i.e., there exists a constant $f_M > 0$ such that

$$0 \le f(x, \xi_2) - f(x, \xi_1) \le f_M(\xi_2 - \xi_1) \quad \text{for a.e. } x \in \Omega \text{ and all } \xi_2 \ge \xi_1.$$

Under these conditions, there exists a unique solution $u \in H_0^1(\Omega)$ of problem (1.1) with (3.1). Moreover, the operator $\mathcal{R}$ is monotone and Lipschitz continuous, as stated in (1.2), with the simple distance metric

$$\mathrm{dist}(u_\ell, u) = \|\nabla(u_\ell - u)\|.$$

The proof of the statements above is postponed to Proposition 7.1, since it uses some well-known arguments and it is not the main focus of the paper.

**Example 3.2 (Systems modeled by (3.1)).** Systems having the flux function $\sigma$ of the form

$$\sigma(x, y) = A(|y|)y, \qquad (3.2)$$

where $A: [0, \infty) \to (0, \infty)$, satisfy (A1) if and only if $0 < \sigma_M \le (A(t))' \le \sigma_M < \infty$ for almost all $t \ge 0$, see, e.g., [33, Appendix A] and the references therein. In this case

$$|A(|y|)y - A(|z|)z| \le \sigma_M |y - z| \quad \text{and} \quad (A(|y|)y - A(|z|)z) \cdot (y - z) \ge \sigma_M |y - z|^2.$$
A common example of this type of system is the mean-curvature flow where

$$A(\varrho) = \sigma_m + (\sigma_M - \sigma_m)/\sqrt{1 + \varrho^2}$$

for the two constants $0 < \sigma_m \le \sigma_M < \infty$. Numerical results are presented for this system in Section 6.1. Examples of other problems close to this category are p-Laplace equations and compressive and non-Newtonian flow equations. However, these satisfy (A1) only under certain additional assumptions on boundedness.

## 3.2 Gradient-independent diffusivity

The second class of problems that we consider foregoes the nonlinear dependence on the gradient. In this case, we assume the following form of $\mathcal{R}$: for $\xi \in H_0^1(\Omega)$,

$$\langle \mathcal{R}(\xi), \varphi \rangle := \tau(\bar{\mathbf{K}}(\cdot))(\mathcal{D}(\cdot, \xi)\nabla\xi + \mathbf{q}(\cdot, \xi)), \nabla\varphi) + \langle f(\cdot, \xi), \varphi \rangle, \quad \varphi \in H_0^1(\Omega). \quad (3.3)$$

We suppose that the auxiliary functions satisfy the following general conditions

**Assumption 3.3 (Gradient-independent diffusivity).**

(B1) The diffusion coefficient $\mathcal{D}: \Omega \times \mathbb{R} \to (0, \infty)$ is bounded and Lipschitz-continuous with respect to all its arguments, i.e., there exists a constant $\mathcal{D}_M > 0$ such that $0 < \mathcal{D} \le \mathcal{D}_M$, and

$$|\mathcal{D}(\mathbf{x}_1, \xi_1) - \mathcal{D}(\mathbf{x}_2, \xi_2)| \le \mathcal{D}_M(|\mathbf{x}_1 - \mathbf{x}_2| + |\xi_1 - \xi_2|) \quad \forall \mathbf{x}_1, \mathbf{x}_2 \in \Omega \text{ and } \xi_1, \xi_2 \in \mathbb{R}.$$

(B2) $f$ satisfies (A2).

(B3) The permeability tensor $\bar{\mathbf{K}}: \Omega \to \mathbb{R}^{d \times d}$ is uniformly symmetric positive definite and bounded, i.e., there exist real constants $K_M \ge K_m > 0$ such that for almost all $\mathbf{x} \in \Omega$,

$$K_M |\mathbf{y}|^2 \le \mathbf{y}^T \bar{\mathbf{K}}(\mathbf{x}) \mathbf{y} \le K_M |\mathbf{y}|^2 \quad \forall \mathbf{y} \in \mathbb{R}^d.$$

(B4) The advective flux $\mathbf{q}: \Omega \times \mathbb{R} \to \mathbb{R}^d$ is bounded by a constant $q_M > 0$, i.e., $|\mathbf{q}| < q_M$. Moreover, it is small with respect the diffusion coefficient $\mathcal{D}$ in that for a constant $\gamma > 0$, the following inequality is satisfied for all $\xi, \zeta \in \mathbb{R}$ and almost all $\mathbf{x} \in \Omega$,

$$\left| \bar{\mathbf{K}}^{\frac{1}{2}}(\mathbf{x}) \left( \mathbf{q}(\mathbf{x}, \zeta) - \mathbf{q}(\mathbf{x}, \xi) - \sum_{j=1}^{d} e_j \int_{\xi}^{\zeta} \partial_{x_j} \mathcal{D}(\mathbf{x}, \varrho) \, d\varrho \right) \right|^2 \le \gamma (f(\mathbf{x}, \zeta) - f(\mathbf{x}, \xi)) \int_{\xi}^{\zeta} \mathcal{D}(\mathbf{x}, \varrho) \, d\varrho.$$

(B5) $\tau$ is a positive real parameter.

The existence and uniqueness of a solution $u \in L^\infty(\Omega) \cap H_0^1(\Omega)$ of the model (1.1) with (3.3) under Assumption 3.3 is shown in Proposition 7.2. These assumptions further ensure that the operator $\mathcal{R}$ is monotone and Lipschitz continuous in the sense of (1.2) with a more complex distance metric

$$\mathrm{dist}(u_\ell, u) = \left\| \bar{\mathbf{K}}^{\frac{1}{2}}(\cdot) \nabla \int_{u}^{u_\ell} \mathcal{D}(\cdot, \varrho) \, d\varrho \right\|.$$

**Example 3.4 (Systems modeled by (3.3)).** We mention two important classes of equations that are modeled by (3.3).

Class I: *Semilinear equations*, corresponding to

$$\mathcal{D} = 1, \quad \bar{\mathbf{K}} = \mathbb{I}, \quad \tau = 1, \text{ and } \mathbf{q}(\mathbf{x}, \xi) = \mathbf{0}. \quad (3.4)$$
(These would of course also fit (3.1) with $\sigma(x, \nabla\xi) = \nabla\xi$). They are used in the modelling of ignition of gases, gravitational influences on stars, in quantum field theory, and in many other applications (see, e.g., [49] and the references therein).

**Class II:** elliptic problems arising from **implicit time-discretization of nonlinear advection-reaction-diffusion equations**, with $\tau > 0$ being the corresponding time-step [52]. Examples are the Fischer-KPP equation, equations used to model flow through porous media (the porous medium equation, the Richards equation), and biofilms [48].

In all cases, Assumption 3.3 is satisfied either directly (semilinear equations), or under some conditions on the auxiliary functions (time-discretized systems). A specific example of the Richards equation is discussed in more detail in Section 6.2.

**Remark 3.5 (Degeneracy).** Conditions (B1)-(B5) allow for degeneracy, i.e., they allow that $D(x, \xi) \searrow 0$ as $\xi \to \pm\infty$, leading to the loss of ellipticity of the problem (3.3), see also Proposition 7.2. This is not covered in the gradient-dependent diffusivity case (3.1).

# 4 Finite element discretization and iterative linearization

To obtain a numerical approximation of $u \in H_0^1(\Omega)$ in (1.1), one usually first discretizes the system in space and then employs an iterative linearization scheme to solve the resulting system of nonlinear algebraic equations. We now describe this procedure in detail.

## 4.1 Finite element discretization

Let $\mathcal{T}_\ell$ denote a triangulation of $\Omega$ with matching and uniformly shape-regular simplicial meshes. Let $\mathcal{V}_\ell$ be the corresponding set of vertices. For each element $K \in \mathcal{T}_\ell$, let $h_K := \text{diam}\{K\}$ denote the diameter of $K$. Let

$$V_{\ell} := H_{0}^{1}(\Omega) \cap \mathcal{P}_{p}(\mathcal{T}_{\ell}), \qquad (4.1)$$

where $\mathcal{P}_p(\mathcal{T}_\ell)$ denotes the space of piecewise polynomials of degree $p \ge 1$ for the triangulation $\mathcal{T}_\ell$. We note that, for all $q \ge 2$,

$$V_{\ell} \subset W_{0}^{1,\infty}(\Omega) \subset L^{\infty}(\Omega) \cap W_{0}^{1,q}(\Omega) \subset L^{\infty}(\Omega) \cap H_{0}^{1}(\Omega), \qquad (4.2)$$

as $V_\ell$ is finite dimensional. Then, $u_\ell \in V_\ell$ is the finite element approximation to the nonlinear problem (1.1) if it satisfies

$$\langle \mathcal{R}(u_\ell), \varphi_\ell \rangle = 0 \quad \forall \varphi_\ell \in V_\ell. \qquad (4.3)$$

To solve the system of nonlinear algebraic equations arising from (4.3), iterative linearization schemes are commonly used, where a sequence $\{u_\ell^i\}_{i \ge 1} \in V_\ell$ is constructed iteratively from an initial guess $v_\ell^0 \in V_\ell$. It yields $u_\ell$ in the limit when the linearization converges. Below, we give an abstract definition of such a linearization scheme.

## 4.2 Iteration-dependent scalar product

For a linearization iteration index $i \ge 0$, let $u_\ell^i \in V_\ell$ be given. We define

$$((\xi, \zeta))_{u_\ell^i} := (L^i \xi, \zeta) + (a^i \nabla \xi, \nabla \zeta), \quad \xi, \zeta \in H_0^1(\Omega) \qquad (4.4)$$

and assume:

**Assumption 4.1 (Coefficients $L^i$ and $a^i$).** The function $L^i: \Omega \to \mathbb{R}$ and the symmetric tensor $a^i: \Omega \to \mathbb{R}^d \times \mathbb{R}^d$ are defined from the data of (3.1) or (3.3) and $u_\ell^i$ and satisfy for functions
0 ≤ L_m^i ≤ L_M^i < ∞ and 0 < a_m^i ≤ a_M^i < ∞, piecewise constant with respect to mesh T_ℓ and sharpest possible, that

$$L_m^i|_K \le L^i(\mathbf{x}) \le L_M^i|_K, \qquad \text{for a.e. } \mathbf{x} \in K, \forall K \in \mathcal{T}_\ell, \qquad (4.5a)$$

$$a_m^i|_K |\mathbf{y}|^2 \le \mathbf{y}^T \mathbf{a}^i(\mathbf{x}) \mathbf{y} \le a_M^i|_K |\mathbf{y}|^2, \quad \forall \mathbf{y} \in \mathbb{R}^d, \text{ for a.e. } \mathbf{x} \in K, \forall K \in \mathcal{T}_\ell. \qquad (4.5b)$$

Consequently,

$$((\xi, \zeta))_{u_{\ell}^{i}} \text{ is a scalar product over } H_{0}^{1}(\Omega) \text{ equivalent to } (\xi, \zeta) + (\nabla\xi, \nabla\zeta) \text{ and } (\nabla\xi, \nabla\zeta). \quad (4.6)$$

**Remark 4.2** (Iteration-dependent scalar product). The scalar product $((\xi, \zeta))_{u_{\ell}^i}$ from (4.4) is of reaction-diffusion type when $L^i \neq 0$ and of pure diffusion type when $L^i = 0$. The first setting in general only appears when the function $f(\mathbf{x}, \xi)$ in (3.1) or (3.3) depends (linearly or nonlinearly) on the second argument $\xi$, i.e., when it represents a (linear or nonlinear) reaction term. When it only depends on the space coordinate $\mathbf{x}$, forming a source term $f(\mathbf{x})$, typically $L^i = 0$, see Tables 1 and 2 below. The scalar product $((\xi, \zeta))_{u_{\ell}^i}$ is always of pure diffusion type for the Zarantonello linearization of Table 1, where $L^i = 0$ by definition.

We will now use this scalar product to define the iterative linearization.

## 4.3 Iterative linearization

Let $u_{\ell}^0 \in V_{\ell}$ be given. For $i \ge 0$, let $u_{\ell}^i \in V_{\ell}$ denote the $i^{th}$ iterate. Then, the $(i+1)^{th}$ iterate $u_{\ell}^{i+1} \in V_{\ell}$ is computed by solving the linear reaction-diffusion or diffusion problem (see Remark 4.2), with the difference of the iterates (increment) in the scalar product (4.4) on the left-hand side and the residual (1.1) on the right-hand side,

$$((u_{\ell}^{i+1} - u_{\ell}^{i}, \varphi_{\ell}))_{u_{\ell}^{i}} = -\langle \mathcal{R}(u_{\ell}^{i}), \varphi_{\ell} \rangle \quad \forall \varphi_{\ell} \in V_{\ell}. \qquad (4.7)$$

Since, recalling (4.2), $u_{\ell}^i \in H_0^1(\Omega)$, problem (4.7) is well posed by the Riesz representation theorem in view of Assumption 4.1 and (4.6). Examples of $((\cdot, \cdot))_{u_{\ell}^i}$ corresponding to the most commonly used linearization schemes are given in Section 4.6, whereas the consistency of the linearization iterations is discussed in a general functional analytic setting in Appendix A.

When expanded using (4.4), problem (4.7) writes: find $u_{\ell}^{i+1} \in V_{\ell}$ such that

$$(L^i u_{\ell}^{i+1}, \varphi_{\ell}) + (\mathbf{a}^i \nabla u_{\ell}^{i+1}, \nabla \varphi_{\ell}) = - \langle \mathcal{R}(u_{\ell}^i), \varphi_{\ell} \rangle + (L^i u_{\ell}^i, \varphi_{\ell}) \\ + (\mathbf{a}^i \nabla u_{\ell}^i, \nabla \varphi_{\ell}) \quad \forall \varphi_{\ell} \in V_{\ell}. \qquad (4.8)$$

Furthermore, we can always identify a scalar-valued source term $S^i \in L^2(\Omega)$ and a vector-valued flux term $\mathcal{F}^i \in L^2(\Omega; \mathbb{R}^d)$, so that (4.8) further rewrites as: find $u_{\ell}^{i+1} \in V_{\ell}$ such that

$$(L^i u_{\ell}^{i+1}, \varphi_{\ell}) + (\mathbf{a}^i \nabla u_{\ell}^{i+1}, \nabla \varphi_{\ell}) = -(S^i, \varphi_{\ell}) - (\mathcal{F}^i, \nabla \varphi_{\ell}) \quad \forall \varphi_{\ell} \in V_{\ell}. \qquad (4.9)$$

Indeed, there holds, for all $\varphi \in H_0^1(\Omega)$,

$$\langle \mathcal{R}(u_{\ell}^{i}), \varphi \rangle - (L^{i} u_{\ell}^{i}, \varphi) - (\mathbf{a}^{i} \nabla u_{\ell}^{i}, \nabla \varphi) = (S^{i}, \varphi) + (\mathcal{F}^{i}, \nabla \varphi), \qquad (4.10)$$

where, in the gradient-dependent case (3.1),

$$S^i := f(\cdot, u_{\ell}^i) - L^i u_{\ell}^i, \qquad (4.11a)$$

$$\mathcal{F}^i := \sigma(\cdot, \nabla u_{\ell}^i) - \mathbf{a}^i \nabla u_{\ell}^i, \qquad (4.11b)$$

whereas in the gradient-independent case (3.3),

$$S^i := f(\cdot, u_{\ell}^i) - L^i u_{\ell}^i, \qquad (4.12a)$$

$$\mathcal{F}^i := \tau \bar{\mathcal{K}}(\cdot) \mathbf{q}(\cdot, u_{\ell}^i) + (\tau \bar{\mathcal{K}}(\cdot) \mathcal{D}(\cdot, u_{\ell}^i) - \mathbf{a}^i) \nabla u_{\ell}^i. \qquad (4.12b)$$
## 4.4 Continuous-level counterpart of (4.7)

Only the discrete problems (4.7) need to be resolved in practice. For our theoretical developments below, though, the following continuous-level counterpart of (4.7) will be useful: we define $u_{\langle \ell \rangle}^{i+1} \in H_0^1(\Omega)$ to be the solution of

$$ ((u_{\langle \ell \rangle}^{i+1} - u_{\ell}^{i}, \varphi))_{u_{\ell}^{i}} = -\langle \mathcal{R}(u_{\ell}^{i}), \varphi \rangle \quad \forall \varphi \in H_0^1(\Omega). \quad (4.13) $$

This problem is again well posed by the Riesz representation theorem. Now $u_{\ell}^{i+1}$ from (4.7) is the finite element approximation of $u_{\langle \ell \rangle}^{i+1}$ from (4.13). As for (4.9), (4.13) can be equivalently written as the linear (reaction-)diffusion problem of finding $u_{\langle \ell \rangle}^{i+1} \in H_0^1(\Omega)$ such that

$$ (L^i u_{\langle \ell \rangle}^{i+1}, \varphi) + (a^i \nabla u_{\langle \ell \rangle}^{i+1}, \nabla \varphi) = -(S^i, \varphi) - (\mathcal{F}^i, \nabla \varphi) \quad \forall \varphi \in H_0^1(\Omega). \quad (4.14) $$

## 4.5 Iteration-dependent norm

In our developments, we crucially rely on the *iteration-dependent norm* associated to the scalar product (4.4) from the iterative linearization (4.7), i.e., given as

$$ \|\xi\|_{1, u_{\ell}^{i}} := ((\xi, \xi))_{u_{\ell}^{i}}^{\frac{1}{2}}, \quad \xi \in H_{0}^{1}(\Omega). \quad (4.15) $$

The Cauchy-Schwarz inequality then reads, for any $\xi, \zeta \in H_0^1(\Omega)$,

$$ ((\xi, \zeta))_{u_{\ell}^{i}} \le \|\xi\|_{1, u_{\ell}^{i}} \|\zeta\|_{1, u_{\ell}^{i}}, \quad (4.16) $$

and we have, for any $\xi \in H_0^1(\Omega)$, the duality expression

$$ \|\xi\|_{1, u_{\ell}^{i}} = \sup_{\varphi \in H_{0}^{1}(\Omega)} \frac{((\xi, \varphi))_{u_{\ell}^{i}}}{\|\varphi\|_{1, u_{\ell}^{i}}}, \quad (4.17) $$

with the equality occurring when $\varphi = \xi/\|\xi\|_{1, u_{\ell}^{i}}$. We will also employ the corresponding dual norm

$$ \|\zeta\|_{-1, u_{\ell}^{i}} := \sup_{\varphi \in H_{0}^{1}(\Omega)} \frac{\langle \zeta, \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^{i}}}, \quad \zeta \in H^{-1}(\Omega). \quad (4.18) $$

Finally, for an open Lipschitz subdomain $\omega \subset \Omega$, we also introduce the local counterparts of $\|\cdot\|_{\pm 1, u_{\ell}^{i}}$ for $\xi \in H_0^1(\Omega)$ and $\zeta \in H^{-1}(\Omega)$ as

$$ \|\xi\|_{1, u_{\ell}^{i}, \omega} := ((L^i \xi, \xi)_\omega + (a^i \nabla \xi, \nabla \xi)_\omega)^{\frac{1}{2}}, \quad (4.19a) $$

$$ \|\zeta\|_{-1, u_{\ell}^{i}, \omega} := \sup_{\varphi \in H_{0}^{1}(\omega)} \frac{\langle \zeta, \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^{i}, \omega}}. \quad (4.19b) $$

## 4.6 Examples of common linearization schemes

We now give examples of some common linearization schemes defining the iteration-dependent scalar product (4.4) and consequently the iterates (4.7).
<table><caption>Table 1: Some common linearizations to solve problem (1.1) with the gradient-dependent expression (3.1) for $\mathcal{R}$. For the Kačanov and Newton iterations, $\sigma(\mathbf{x}, \mathbf{y}) = A(|\mathbf{y}|)\mathbf{y}$ is assumed, as in Example 3.2.</caption><thead><tr><th>Scheme</th><th>L<sup>i</sup>(<b>x</b>)</th><th>a<sup>i</sup>(<b>x</b>)</th></tr></thead><tbody><tr><td>Kačanov [68, Chapter 25]</td><td>&part;<sub>&xi;</sub>f(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td><td>A(|&nabla;u<sub>&ell;</sub><sup>i</sup>|)</td></tr><tr><td>Newton [22]</td><td>&part;<sub>&xi;</sub>f(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td><td>A(|&nabla;u<sub>&ell;</sub><sup>i</sup>|) + <sup>A'(|&nabla;u<sub>&ell;</sub><sup>i</sup>|)</sup>&frasl;<sub>|&nabla;u<sub>&ell;</sub><sup>i</sup>|</sub> &nabla;u<sub>&ell;</sub><sup>i</sup> &otimes; &nabla;u<sub>&ell;</sub><sup>i</sup></td></tr><tr><td>Zarantonello [66]</td><td>0</td><td>&Lambda; (constant) &gt; 0</td></tr></tbody></table>

### 4.6.1 Gradient-dependent diffusivity

For the gradient-dependent diffusivity case (3.1), the three most standard linearization schemes are presented in Table 1. The Kačanov (fixed point) iterations are considered in Chapter 25 of [68], and the convergence of the scheme is proven under certain restrictions on $\sigma(\mathbf{x}, \mathbf{y})$. The Newton [22] linearization adds a first-order correction and converges quadratically for a sufficiently precise initial guess. The Zarantonello [66] scheme converges linearly irrespective of the initial guess, for a constant $\Lambda \ge \sigma_M^2/\sigma_m$, where $\sigma_M$ and $\sigma_m$ stem from (A1). We refer to [33,36] and the references therein for additional properties.

From (4.11) and Table 1, one in particular obtains that the linearizations (4.7) read here as (4.9) with

$$ S^i = f(\cdot, u_{\ell}^i) - \partial_{\xi} f(\cdot, u_{\ell}^i) u_{\ell}^i \quad \text{and} \quad \mathcal{F}^i = 0 \qquad [\text{Kačanov}], \quad (4.20a) $$

$$ S^i = f(\cdot, u_{\ell}^i) - \partial_{\xi} f(\cdot, u_{\ell}^i) u_{\ell}^i \quad \text{and} \quad \mathcal{F}^i = -A'(|\nabla u_{\ell}^i|) |\nabla u_{\ell}^i| |\nabla u_{\ell}^i| \quad [\text{Newton}], \quad (4.20b) $$

$$ S^i = f(\cdot, u_{\ell}^i) \quad \text{and} \quad \mathcal{F}^i = \sigma(\cdot, \nabla u_{\ell}^i) - \Lambda \nabla u_{\ell}^i \qquad [\text{Zarantonello}]. \quad (4.20c) $$

### 4.6.2 Gradient-independent diffusivity

For the gradient-independent diffusivity case (3.3), the commonly used linearization schemes are collected in Table 2. The four schemes presented are all linearly converging under some restrictions on the parameter $\tau > 0$.

<table><caption>Table 2: Some common linearization techniques to solve the problem (1.1) with the gradient-independent expression (3.3) for $\mathcal{R}$.</caption><thead><tr><th>Scheme</th><th>L<sup>i</sup>(<b>x</b>)</th><th>a<sup>i</sup>(<b>x</b>)</th></tr></thead><tbody><tr><td>Picard [13]</td><td>&part;<sub>&xi;</sub>f(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td><td>&tau;K&#773;(<b>x</b>) D(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td></tr><tr><td>Jäger-Kačur [39]</td><td>max<sub>&xi; &isin; &#8477;</sub> (<sup>f(<b>x</b>,&xi;)-f(<b>x</b>,u<sub>&ell;</sub><sup>i</sup>)</sup>&frasl;<sub>&xi;-u<sub>&ell;</sub><sup>i</sup></sub>)</td><td>&tau;K&#773;(<b>x</b>) D(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td></tr><tr><td>L-scheme [50, 55]</td><td>L (constant) &ge; <sup>1</sup>&frasl;<sub>2</sub> sup &part;<sub>&xi;</sub>f</td><td>&tau;K&#773;(<b>x</b>) D(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td></tr><tr><td>M-scheme [52]</td><td>&part;<sub>&xi;</sub>f(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>) + M&tau; (constant)</td><td>&tau;K&#773;(<b>x</b>) D(<b>x</b>, u<sub>&ell;</sub><sup>i</sup>)</td></tr></tbody></table>

**Remark 4.3** (Newton(–Raphson) linearization). For the case (3.3), the Newton(–Raphson) iterative scheme does not yield a symmetric linearization problem in general. We will discuss it separately in Appendix C.

From (4.12) and Table 2, one in particular obtains that the linearizations (4.7) read here as (4.9) with

$$ S^i = f(\cdot, u_{\ell}^i) - L^i u_{\ell}^i \quad \text{and} \quad \mathcal{F}^i = \tau \bar{K}(\cdot) q(\cdot, u_{\ell}^i). \qquad (4.21) $$

Observe that, for all the schemes in Table 2, the second term cancels out in (4.12b).
# 5 Main results

In this section, we show that there exists an orthogonality relation between the total error in the finite element solution $u_{\ell}^i$, the linearization error, and the discretization error. The iteration-dependent norm defined in Section 4.5 is crucial for this. We further state how this result allows us to estimate the total error using *guaranteed* (binding above), *efficient* (binding below, up to a constant), and *fully computable* a posteriori estimators *robust* with respect to the strength of the nonlinearities in (3.1) or (3.3). Congruently, the distinction of the error components allows us to balance them and to stop the iterative linearizations timely.

## 5.1 An orthogonal decomposition of the error

We start by the orthogonal decomposition of the total error in the finite element approximation $u_{\ell}^i$.

### 5.1.1 Discretization residual and discretization error

Stemming from the linearized problems (4.7) and (4.13), we define the discretization residual $\mathcal{R}_{\text{disc}}^{u_{\ell}^i}(\xi) \in H^{-1}(\Omega)$, for $\xi \in H_0^1(\Omega)$, as

$$ \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(\xi), \varphi \rangle := \langle (\xi - u_{\ell}^i, \varphi) \rangle_{u_{\ell}^i} + \langle \mathcal{R}(u_{\ell}^i), \varphi \rangle, \quad \varphi \in H_0^1(\Omega). \quad (5.1) $$

In particular, using the definition (4.4) and (4.10), with the notations (4.11)–(4.12),

$$ \begin{aligned} \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle &= \langle (u_{\ell}^{i+1} - u_{\ell}^i, \varphi) \rangle_{u_{\ell}^i} + \langle \mathcal{R}(u_{\ell}^i), \varphi \rangle \\ &= (L^i u_{\ell}^{i+1}, \varphi) + (a^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i, \nabla \varphi) + (S^i, \varphi). \end{aligned} \quad (5.2) $$

Note that, from (4.13), there holds

$$ \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\langle \ell \rangle}^{i+1}), \varphi \rangle = 0 \quad \forall \varphi \in H_0^1(\Omega), \quad (5.3) $$

which justifies the name “residual”: when the finite-dimensional space $V_{\ell}$ approaches $H_0^1(\Omega)$, $\mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1})$ vanishes. More precisely, since $u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1}$ is the difference between the discrete and the continuous level solutions from (4.7)–(4.13), using the following result, we find that $\|\mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1})\|_{-1, u_{\ell}^i}$ represents the discretization error of the linearization step $i+1$:

**Lemma 5.1** (Discretization residual–discretization error relation). There holds

$$ \| \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}) \|_{-1, u_{\ell}^i} = \| u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1} \|_{1, u_{\ell}^i}. $$

*Proof.* From (5.1) and (5.3), we have that

$$ \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle \stackrel{(5.3)}{=} \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}) - \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\langle \ell \rangle}^{i+1}), \varphi \rangle \stackrel{(5.1)}{=} \langle (u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1}, \varphi) \rangle_{u_{\ell}^i} \quad \forall \varphi \in H_0^1(\Omega). $$

Thus, the result follows from (4.17)–(4.18). □
### 5.1.2 Total residual and total error

In order to characterize the total error in $u_{\ell}^{i} \in V_{\ell}$, we need to consider the total residual $\mathcal{R}(u_{\ell}^{i})$ stemming from (1.1); as motivated in Section 1, we consider $\|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}$ as the total error measure. Observe from (4.18) that if $u \in H_{0}^{1}(\Omega)$ solves (1.1), then

$$ \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} = \sup_{\varphi \in H_{0}^{1}(\Omega)} \frac{\langle \mathcal{R}(u_{\ell}^{i}), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^{i}}} = \sup_{\varphi \in H_{0}^{1}(\Omega)} \frac{\langle \mathcal{R}(u_{\ell}^{i}) - \mathcal{R}(u), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^{i}}}, \quad (5.4) $$

which is, due to (1.2) and (4.6), equivalent to the error metric $\mathrm{dist}(u_{\ell}^{i}, u)$. Additionally, using

$$ \langle \mathcal{R}(u_{\ell}^{i}), \varphi \rangle \stackrel{(4.13)}{=} \left( \|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}, \varphi\| \right)_{u_{\ell}^{i}} \quad \forall \varphi \in H_{0}^{1}(\Omega), \quad (5.5) $$

and employing (4.17)–(4.18), one has:

**Lemma 5.2 (Total residual–total error relation).** There holds

$$ \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} = \|\|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}\|\|_{1, u_{\ell}^{i}}. $$

### 5.1.3 An orthogonal decomposition of the total residual/error

The remaining component of the total residual/error is the linearization error; it comes naturally out of the following decomposition result:

**Theorem 5.3 (Orthogonal decomposition of the total residual/error into linearization and discretization components).** Let the nonlinear elliptic operator $\mathcal{R}: H_{0}^{1}(\Omega) \to H^{-1}(\Omega)$ be introduced in (1.1) and satisfy (1.2). Let the finite element linearization sequence $\{u_{\ell}^{i}\}_{i \ge 0} \subset V_{\ell}$ be given by (4.7), for the iteration-dependent scalar product (4.4) stemming from the given linearization and satisfying Assumption 4.1. Let the affine elliptic operator $\mathcal{R}_{\mathrm{disc}}: H_{0}^{1}(\Omega) \to H^{-1}(\Omega)$ be given by (5.1), stemming from (4.13). Let finally the norms $\|\cdot\|_{\pm 1, u_{\ell}^{i}}$ be given in (4.15), (4.18). Then, for all linearization steps $i \ge 0$,

$$ \underbrace{\|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2}}_{\text{total residual/error}} = \underbrace{\|u_{\ell}^{i} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}}^{2}}_{\text{linearization error}} + \underbrace{\|\mathcal{R}_{\mathrm{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})\|_{-1, u_{\ell}^{i}}^{2}}_{\substack{\text{discretization error} \\ \text{residual}}} \cdot \|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}\|_{1, u_{\ell}^{i}} $$

*Proof.* One has

$$ \begin{align*} \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2} & \stackrel{\text{Lemma 5.2}}{=} \|\|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}\|\|_{1, u_{\ell}^{i}}^{2} = \|\|(u_{\ell}^{i} - u_{\ell}^{i+1}) + (u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1})\|\|_{1, u_{\ell}^{i}}^{2} \\ & \stackrel{(4.15)}{=} \|\|u_{\ell}^{i} - u_{\ell}^{i+1}\|\|_{1, u_{\ell}^{i}}^{2} + \|\|u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1}\|\|_{1, u_{\ell}^{i}}^{2} + 2\|(u_{\ell}^{i} - u_{\ell}^{i+1}, u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1})\|_{u_{\ell}^{i}} \\ & \stackrel{(4.7),(4.13)}{=} \|\|u_{\ell}^{i} - u_{\ell}^{i+1}\|\|_{1, u_{\ell}^{i}}^{2} + \|\|u_{\ell}^{i+1} - u_{\langle \ell \rangle}^{i+1}\|\|_{1, u_{\ell}^{i}}^{2} \\ & \stackrel{\text{Lemma 5.1}}{=} \|\|u_{\ell}^{i} - u_{\ell}^{i+1}\|\|_{1, u_{\ell}^{i}}^{2} + \|\|\mathcal{R}_{\mathrm{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})\|\|_{-1, u_{\ell}^{i}}^{2}, \end{align*} $$

where the penultimate equality follows from the Galerkin orthogonality between $u_{\ell}^{i+1}$ and $u_{\langle \ell \rangle}^{i+1}$, using that $u_{\ell}^{i} - u_{\ell}^{i+1}$ lies in the finite element space $V_{\ell}$. $\square$
## 5.2 Guaranteed and locally efficient a posteriori error bounds robust with respect to the strength of the nonlinearity

Based on Theorem 5.3, estimating $||\mathcal{R}(u_{\ell}^i)||_{-1, u_{\ell}^i}$ amounts to estimating $||\mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1})||_{-1, u_{\ell}^i}$, where the latter corresponds to the numerical approximation (4.7) of the linear (reaction-)diffusion problem (4.13). For this purpose, we will follow [1, 2, 42, 43, 60, 64], the pure diffusion approaches [10, 20, 21, 24], but most precisely [60, Definition 2.3].

For each vertex $\mathbf{a}$ in the set of vertices $\mathcal{V}_\ell$, let the patch $\mathcal{T}_\mathbf{a}$ be the set of mesh elements from $\mathcal{T}_\ell$ sharing $\mathbf{a}$. We denote by $\psi_\mathbf{a}$ the corresponding hat function (piecewise affine function with value 1 in $\mathbf{a}$ and 0 in the other vertices), with the subdomain $\omega_\mathbf{a}$ being the interior of the support of $\psi_\mathbf{a}$. Let $h_{\omega_\mathbf{a}}$ be the diameter of $\omega_\mathbf{a}$.

As in (4.1), let $\mathcal{P}_p(K)$ stand for the space of polynomials of total degree at most $p$ on the simplex $K \in \mathcal{T}_\ell; \mathcal{P}_p(K; \mathbb{R}^d)$ denotes its componentwise extension to $\mathbb{R}^d$. The broken (elementwise) spaces $\mathcal{P}_p(\mathcal{T}_\mathbf{a})$ and $\mathcal{RT}_p(\mathcal{T}_\mathbf{a})$ on the patch $\mathcal{T}_\mathbf{a}$ are then defined by, cf. [7],

$$ \mathcal{P}_p(\mathcal{T}_\mathbf{a}) := \{v_\ell \in L^2(\omega_\mathbf{a}), \quad v_\ell|_K \in \mathcal{P}_p(K) \quad \forall K \in \mathcal{T}_\mathbf{a}\}, \\ \mathcal{RT}_p(\mathcal{T}_\mathbf{a}) := \{v_\ell \in L^2(\omega_\mathbf{a}; \mathbb{R}^d), \quad v_\ell|_K \in \mathcal{RT}_p(K) := \mathcal{P}_p(K; \mathbb{R}^d) + x\mathcal{P}_p(K) \quad \forall K \in \mathcal{T}_\mathbf{a}\}. $$

We remark that $\mathcal{P}_p(\mathcal{T}_\mathbf{a}; \mathbb{R}^d) \subsetneq \mathcal{RT}_p(\mathcal{T}_\mathbf{a}) \subsetneq \mathcal{P}_{p+1}(\mathcal{T}_\mathbf{a}; \mathbb{R}^d)$. Decompose $\mathcal{V}_\ell$ into interior vertices $\mathcal{V}_\ell^{\text{int}}$ and boundary vertices $\mathcal{V}_\ell^{\text{ext}}$. The local mixed finite element spaces are then defined by

$$ Q_{\mathbf{a}} := \begin{cases} \mathcal{P}_p(\mathcal{T}_{\mathbf{a}}) & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{ext}} \text{ or } L_{\text{M}}^{i}|_{\omega_{\mathbf{a}}} > 0, \\ \{v_{\ell} \in \mathcal{P}_p(\mathcal{T}_{\mathbf{a}}), \quad (v_{\ell}, 1)_{\omega_{\mathbf{a}}} = 0\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}} \text{ and } L_{\text{M}}^{i}|_{\omega_{\mathbf{a}}} = 0, \end{cases} \quad (5.6a) $$

$$ V_{\mathbf{a}} := \begin{cases} \{v_{\ell} \in \mathcal{RT}_p(\mathcal{T}_{\mathbf{a}}) \cap \mathcal{H}(\text{div}, \omega_{\mathbf{a}}), \quad v_{\ell} \cdot n = 0 \text{ on } \partial \omega_{\mathbf{a}}\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}}, \\ \{v_{\ell} \in \mathcal{RT}_p(\mathcal{T}_{\mathbf{a}}) \cap \mathcal{H}(\text{div}, \omega_{\mathbf{a}}), \quad v_{\ell} \cdot n = 0 \text{ on } \partial \omega_{\mathbf{a}} \setminus \{\psi_{\mathbf{a}} > 0\}\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{ext}}. \end{cases} \quad (5.6b) $$

In order to define our estimators, we solve independent small discrete local linear reaction-diffusion problems in a dual formulation on the vertex patches $\mathcal{T}_\mathbf{a}$. They ultimately yield the equilibrated flux-potential pair $\sigma_\ell^{i+1}$, $\phi_\ell^{i+1}$ and read (details of the notation are postponed to Section 8.1 below):

**Definition 5.4** (Equilibrated flux $\sigma_\ell^{i+1}$ and potential $\phi_\ell^{i+1}$). Let $u_\ell^{i+1} \in V_\ell$ be defined by (4.7) or, equivalently, (4.9). For each vertex $\mathbf{a} \in \mathcal{V}_\ell$, let $(\sigma_\mathbf{a}^{i+1}, \phi_\mathbf{a}^{i+1}) \in \mathbf{V}_\mathbf{a} \times Q_\mathbf{a}$ be defined by the constrained minimization problem

$$ (\sigma_{\mathbf{a}}^{i+1}, \phi_{\mathbf{a}}^{i+1}) := \arg \min_{(v_{\ell}, q_{\ell}) \in \mathbf{V}_{\mathbf{a}} \times Q_{\mathbf{a}}, \nabla \cdot \mathbf{v}_{\ell} + \Pi_{\ell,p}(L^i q_{\ell}) = -\nabla \psi_{\mathbf{a}} \cdot \Pi_{\ell,p-1}^{\text{RT}}(a^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i) - \Pi_{\ell,p}(\psi_{\mathbf{a}} S^i)} \left[ \begin{aligned} & (\|w^i(a_{\text{m}}^i)\|_{\omega_{\mathbf{a}}}^{-1} (\psi_{\mathbf{a}} \Pi_{\ell,p-1}^{\text{RT}}(a^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i) + v_\ell)\|^2_{\omega_{\mathbf{a}}} \\ & + \|(L^i)^{\frac{1}{2}}(\Pi_{\ell,p}(\psi_{\mathbf{a}} u_{\ell}^{i+1}) - q_\ell)\|^2_{\omega_{\mathbf{a}}} \end{aligned} \right], \quad (5.7) $$

with, respectively, the coefficients $0 \le L_{\text{M}}^i$, $0 < a_{\text{m}}^i$ defined in (4.5), the scalar- and vector-valued source terms $S^i$ and $\mathcal{F}^i$ in (4.11) or (4.12), the elementwise constant weight function $w^i > 0$ in (8.3a), and the projection operators $\Pi_{\ell,p} : L^2(\Omega) \to \mathcal{P}_p(\mathcal{T}_\ell)$, $\Pi_{\ell,p-1}^{\text{RT}} : L^2(\Omega; \mathbb{R}^d) \to \mathcal{RT}_{p-1}(\mathcal{T}_\ell; \mathbb{R}^d)$ in (8.1) below. Then, after extending $\sigma_\mathbf{a}^{i+1}$, $\phi_\mathbf{a}^{i+1}$ by zero from each patch subdominant $\omega_\mathbf{a}$ to $\Omega$, we define

$$ \sigma_{\ell}^{i+1} := \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \sigma_{\mathbf{a}}^{i+1}, \quad \phi_{\ell}^{i+1} := \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \phi_{\mathbf{a}}^{i+1}. \quad (5.8) $$

**Remark 5.5** (Flux equilibration). The flux equilibration of Definition 5.4 is in general of reaction-diffusion type, following [60, Definition 2.3]. If, however, $L^i = 0$, as in Remark 4.2, Definition 5.4 boils down to the usual pure diffusion equilibration of [10, 20, 21, 24] where the weight function $w^i$ equals one and the potentials $\phi_\mathbf{a}^{i+1}$ and $\phi_\ell^{i+1}$ are not needed.
Definition 5.4 in particular implies:

**Lemma 5.6** (Properties of $\sigma_{\ell}^{i+1}$ and $\phi_{\ell}^{i+1}$). The equilibrated flux and potential from Definition 5.4 are well defined and satisfy

$$ \sigma_{\ell}^{i+1} \in \mathcal{RT}_p(\mathcal{T}_{\ell}) \cap \mathcal{H}(\text{div}, \Omega), \quad \phi_{\ell}^{i+1} \in \mathcal{P}_p(\mathcal{T}_{\ell}). \qquad (5.9) $$

The potential contributions satisfy the mean value property

$$ (L^i \phi_{\mathbf{a}}^{i+1}, 1)_{\omega_{\mathbf{a}}} = (L^i u_{\ell}^{i+1}, \psi_{\mathbf{a}})_{\omega_{\mathbf{a}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}} \qquad (5.10) $$

and the flux fulfils the mass balance

$$ \Pi_{\ell,p}(L^i \phi_{\ell}^{i+1}) + \nabla \cdot \sigma_{\ell}^{i+1} = -\Pi_{\ell,p} S^i. \qquad (5.11) $$

Let a subdomain $\omega \subseteq \Omega$ correspond to some mesh elements from $\mathcal{T}_{\ell}$ and let $\tilde{\omega}$ correspond to $\omega$ and the mesh elements neighboring by vertex. With the equilibrated flux $\sigma_{\ell}^{i+1}$ and potential $\phi_{\ell}^{i+1}$, we can now define our (discretization error) estimators:

$$ \eta_{F,\omega}^i := \|\mathbf{w}^i(a_m^i)^{-\frac{1}{2}}(\Pi_{\ell,p-1}^{\text{RT}}(\mathbf{a}^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i) + \sigma_{\ell}^{i+1})\|_{\omega}, \qquad (5.12a) $$

$$ \eta_{S,\omega}^i := \|(L^i)^{\frac{1}{2}}(u_{\ell}^{i+1} - \phi_{\ell}^{i+1})\|_{\omega}, \qquad (5.12b) $$

$$ \eta_{\text{osc},\omega}^i := \left( \sum_{K \in \mathcal{T}_{\ell}, K \subseteq \omega} (\tilde{w}_K^i \|(I - \Pi_{\ell,p})(S^i + L^i \phi_{\ell}^{i+1})\|_K)^2 \right)^{\frac{1}{2}}, \qquad (5.12c) $$

$$ \eta_{\text{quad,F},\omega}^i := \|(\mathbf{a}^i)^{-\frac{1}{2}}((\mathbb{I} - \Pi_{\ell,p-1}^{\text{RT}})(\mathbf{a}^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i))\|_{\omega}, \qquad (5.12d) $$

$$ \eta_{\text{quad,S},\omega}^i := \left( \sum_{\mathbf{a} \in \mathcal{V}_{\ell}, \omega_{\mathbf{a}} \subseteq \tilde{\omega}} \|(L^i)^{\frac{1}{2}}(\Pi_{\ell,p}(\psi_{\mathbf{a}} u_{\ell}^{i+1}) - \Pi_{\ell,p,L^i}(\psi_{\mathbf{a}} u_{\ell}^{i+1}))\|_{\omega_{\mathbf{a}}} \right)^{\frac{1}{2}}, \qquad (5.12e) $$

$$ \eta_{\text{quad,S,osc},\omega}^i := \left( \sum_{\mathbf{a} \in \mathcal{V}_{\ell}, \omega_{\mathbf{a}} \subseteq \tilde{\omega}} \frac{h_{\omega_{\mathbf{a}}}^2}{\min_{\omega_{\mathbf{a}}} a_m^i} \|(I - \Pi_{\ell,p})(\psi_{\mathbf{a}}(S^i + L^i u_{\ell}^{i+1}))\|_{\omega_{\mathbf{a}}}^2 \right)^{\frac{1}{2}}. \qquad (5.12f) $$

They respectively express the $\mathcal{H}(\text{div}, \Omega)$ non-conformity of the direct finite element fluxes $\mathbf{a}^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i$, a possible weighted discrepancy between $u_{\ell}^{i+1}$ and $\phi_{\ell}^{i+1}$ important for large values of $L^i$, the possible oscillations in the source term $S^i$ and $L^i \phi_{\ell}^{i+1}$, and quadrature errors for the flux and source terms. Note that the quadrature estimator $\eta_{\text{quad,S},\omega}^i$ vanishes if $L^i$ is elementwise constant; similarly, if $\mathbf{a}^i$ is elementwise constant and $\mathcal{F}^i \in \mathcal{P}_p(\mathcal{T}_{\ell}, \mathbb{R}^d)$ (often satisfied in the lowest polynomial degree setting $p=1$), then $\eta_{\text{quad,F},\omega}^i = 0$.

Introducing

$$ \eta_{\text{disc},\omega}^i := \left( \sum_{\substack{K \in \mathcal{T}_{\ell} \\ K \subseteq \omega}} (\eta_{F,K}^i + \eta_{S,K}^i + \eta_{\text{osc},K}^i + \eta_{\text{quad,F,K}}^i)^2 \right)^{\frac{1}{2}}, \qquad (\text{discretization estimator}) \qquad (5.13a) $$

$$ \eta_{\text{lin},\omega}^i := \|u_{\ell}^{i+1} - u_{\ell}^i\|_{1, u_{\ell}^i, \omega}, \qquad (\text{linearization estimator}) \qquad (5.13b) $$

$$ \eta_{\omega}^i := ([\eta_{\text{lin},\omega}^i]^2 + [\eta_{\text{disc},\omega}^i]^2)^{\frac{1}{2}}, \qquad (\text{total estimator}) \qquad (5.13c) $$

our second main result reads:
**Theorem 5.7** (Guaranteed, locally efficient, and robust a posteriori estimates). Let the assumptions of Theorem 5.3 be satisfied. For a subdomain $\omega \subseteq \Omega$ corresponding to some mesh elements from $\mathcal{T}_\ell$, let $\tilde{\omega}$ correspond to $\omega$ and the mesh elements neighboring by vertex. For any linearization step $i \ge 0$, there holds

$$ \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2} \leq [\eta_{\Omega}^{i}]^{2}, \qquad (5.14a) $$
(guaranteed reliability)

$$ [\eta_{\omega}^{i}]^{2} \lesssim \vartheta_{\tilde{\omega}, \ell, i}^{2} \left\| \mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1}) \right\|_{-1, u_{\ell}^{i}, \tilde{\omega}}^{2} + [\eta_{\text{lin}, \omega}^{i}]^{2} + [\eta_{\text{osc}, \tilde{\omega}}^{i}]^{2} + \vartheta_{\tilde{\omega}, \ell, i}^{2} [\eta_{\text{quad}, \text{F}, \tilde{\omega}}^{i}]^{2} + [\eta_{\text{quad}, \text{S}, \tilde{\omega}}^{i}]^{2} \\ + [\eta_{\text{quad}, \text{S}, \text{osc}, \tilde{\omega}}^{i}]^{2}, \qquad (5.14b) $$
(local efficiency)

$$ [\eta_{\Omega}^{i}]^{2} \lesssim \vartheta_{\Omega, \ell, i}^{2} \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2} + [\eta_{\text{osc}, \Omega}^{i}]^{2} + \vartheta_{\Omega, \ell, i}^{2} [\eta_{\text{quad}, \text{F}, \Omega}^{i}]^{2} + [\eta_{\text{quad}, \text{S}, \Omega}^{i}]^{2} + [\eta_{\text{quad}, \text{S}, \text{osc}, \Omega}^{i}]^{2}, \qquad (5.14c) $$
(global efficiency)

where the local variability constant $\vartheta_{\omega, \ell, i} \ge 1$ is defined as $\vartheta_{\omega, \ell, i} := \max_{\mathbf{a} \in \mathcal{V}_\ell, \omega_\mathbf{a} \subseteq \omega} \vartheta_{\mathbf{a}, \ell, i}$ with

$$ \vartheta_{\mathbf{a}, \ell, i}^{2} := \frac{\max_{\omega_{\mathbf{a}}} a_{\mathbf{M}}^{i}}{\min_{\omega_{\mathbf{a}}} a_{\mathbf{m}}^{i}} \qquad \text{if } \frac{h_{\omega_{\mathbf{a}}}^{2} \max_{\omega_{\mathbf{a}}} L_{\mathbf{M}}^{i}}{\max_{\omega_{\mathbf{a}}} a_{\mathbf{M}}^{i}} \le 1, \qquad (5.15a) $$
(diffusion-dominated case)

$$ \vartheta_{\mathbf{a}, \ell, i}^{2} \text{ only depends on } \frac{\max_{\omega_{\mathbf{a}}} a_{\mathbf{M}}^{i}}{\min_{\omega_{\mathbf{a}}} a_{\mathbf{m}}^{i}}, \frac{\max_{\omega_{\mathbf{a}}} L_{\mathbf{M}}^{i}}{\min_{\omega_{\mathbf{a}}} L_{\mathbf{m}}^{i}} \qquad \text{if } \frac{h_{\omega_{\mathbf{a}}}^{2} \max_{\omega_{\mathbf{a}}} L_{\mathbf{M}}^{i}}{\max_{\omega_{\mathbf{a}}} a_{\mathbf{M}}^{i}} > 1 \qquad (5.15b) $$
(reaction-dominated case)

and measures the maximum local patch-wise variation of coefficients $L^i$ and $a^i$ of (4.4) in $\omega$, and the constant in $\lesssim$ only depends on the space dimension $d$, the shape-regularity constant of $\mathcal{T}_\ell$, and the polynomial degree $p$.

The proof of Theorem 5.3 is postponed to Section 8. We conclude this section by three remarks.

**Remark 5.8** (Theorem 5.7). The estimate (5.14a) gives a guaranteed and fully computable upper bound on the error. The error lower bounds (5.14b)-(5.14c) feature a generic constant hidden in $\lesssim$ (only depending on $p$, $d$, and the regularity of $\mathcal{T}_\ell$), so that their quality (robustness) solely hinges on $\vartheta_{\mathbf{a}, \ell, i}$ from (5.15).

**Remark 5.9** (Robustness in the Zarantonello case). The Zarantonello linearization of Table 1 gives $L^i = 0$ and $a^i = \Lambda = \text{constant}$, so that case (5.15a) always applies and the variability constants $\vartheta_{\mathbf{a}, \ell, i} = 1$. In this case, the estimates of Theorem 5.3 are always robust with respect to the strength of the nonlinearity.

**Remark 5.10** (Cases (5.15a) and (5.15b)). Case (5.15a) corresponds to a situation where the reaction coefficients $L^i$ from the linearization scalar product (4.4) do not dominate over the diffusion coefficients $a^i$ from (4.4); note the presence of the mesh size $h_{\omega_\mathbf{a}}$ implying that this case always prevails under sufficient local mesh refinement. Also recall from Remark 4.2 and Tables 1 and 2 that typically $L^i = 0$ when $f(\mathbf{x}, \xi)$ in (3.1) or (3.3) only depends on the space coordinate $\mathbf{x}$, i.e., when $f(\mathbf{x})$ is a source term.

**Remark 5.11** (The variability constants $\vartheta_{\mathbf{a}, \ell, i}$ from (5.15)). The constants $\vartheta_{\mathbf{a}, \ell, i}$ from (5.15) describe a local variation of the linearization coefficient $a^i$ via the piecewise constants $a_{\mathbf{m}}^i, a_{\mathbf{M}}^i$ from Assumption 4.1 in (5.15a), and additionally on the local variation of $L^i$ via $L_{\mathbf{m}}^i, L_{\mathbf{M}}^i$ in (5.15b). This is expected to be (much) smaller than the global variation of $a^i$ and $L^i$ over
the whole computational domain $\Omega$ due to the interior regularity of solutions of elliptic equations [26, Section 6.3], and definitely much smaller than the ratio $\lambda_M/\lambda_m$ from (1.4) or (1.5); $\vartheta_{\alpha,\ell,i}$ can be small and/or tend to one with local mesh refinement even for strong nonlinearities. Moreover, even in the worst case, at least for the gradient-dependent diffusivity (3.1) in the case of Example 3.2 and with f only depending on the space coordinate $x$, $\vartheta_{\alpha,\ell,i} \le (\sigma_M/\sigma_m)^{1/2}$ for the Kačanov and Newton linearization schemes of Table 1 as discussed in [33, Sections 2.3.2 and 3.7], which still outperforms (1.4) and (1.5) by the factor $(\sigma_M/\sigma_m)^{1/2}$.

**Remark 5.12** (A posteriori assessment of robustness). Whenever $\vartheta_{\Omega,\ell,i}$ is of order of unity, the estimates of Theorem 5.7 are robust with respect to the strength of the nonlinearity. Crucially, $\vartheta_{\Omega,\ell,i}$ can be easily computed at any mesh $\mathcal{T}_\ell$ and linearization iteration $i$ in any setting (problem, linearization scheme) in (5.15a), so that robustness is assessed a posteriori.

## 5.3 Balancing of the linearization and discretization components and stopping criteria for iterative linearizations

The orthogonal error decomposition of Theorem 5.3 together with the computable estimators of Theorem 5.7 yield the following adaptive stopping criterion for iterative linearizations that ensures a natural balance between discretization and linearization components, in that it requests that the linearization component only forms a user-given $\mu$-fraction of the total error.

**Algorithm 5.13** (Adaptive linearization). Let $\mu \in (0,1)$ and $u_\ell^0 \in V_\ell$ be fixed. Pursue the iterative linearization (4.7) until, for some $\bar{i} \ge 0$,

$$ \eta_{\text{lin},\Omega}^{\bar{i}} \le \mu \eta_{\Omega}^{\bar{i}}, \qquad (5.16) $$

where $\eta_{\text{lin},\Omega}^{\bar{i}}$ and $\eta_{\Omega}^{\bar{i}}$ are computed following (5.13).

Moreover, the linearization estimator $\eta_{\text{lin},\Omega}^i$ can also be used to switch between linearization schemes in an efficient way, see [61] for an example. However, this is not pursued in this work.

# 6 Numerical experiments

In this section, we test the numerical performance of the estimators of Theorem 5.7 and the stopping criteria of Algorithm 5.13. In what follows, a fixed parameter $\mu = 0.05$ from (5.16) has been used. Three levels of discretization are considered with an initial mesh $\mathcal{T}_1$ fixed and its uniform refinements $\mathcal{T}_\ell$, where we choose $\ell \in \{2, 4\}$. Following Section 4.1, the corresponding finite element spaces $V_\ell$ are given as $\mathcal{P}_1(\mathcal{T}_\ell) \cap H^1(\Omega)$, i.e., $p = 1$. For linearization step $i \ge 0$ and $\ell \in \{1, 2, 4\}$, let $u_{\text{ref}}^{i+1} \in V_{10}$ (i.e. $\ell = 10$) be a fine-mesh approximation to $u_{\ell}^{i+1}$ from (4.13) that we use for approximate evaluation of the error; indeed, by Lemma 5.2, one has that

$$ \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} = \|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}\|_{1, u_{\ell}^{i}} \approx \|u_{\ell}^{i} - u_{\text{ref}}^{i+1}\|_{1, u_{\ell}^{i}}. \qquad (6.1) $$

With these definitions, we express the quality of the estimates by introducing the global and local effectivity indices:

$$ \text{Eff. Ind.} := \eta_{\Omega}^{i} / \|u_{\ell}^{i} - u_{\text{ref}}^{i+1}\|_{1, u_{\ell}^{i}}, \qquad (6.2a) $$

$$ (\text{Eff. Ind.})_K := \eta_K^i / \|u_\ell^i - u_{\text{ref}}^{i+1}\|_{1, u_\ell^i, K}, \quad K \in \mathcal{T}_\ell, \qquad (6.2b) $$

where $\eta_{\Omega}^i$ and $\eta_K^i$ are defined in (5.13c). The numerical solutions are computed using the FreeFem++ software [34] and the codes are available online in a Git repository. We consider the following two test cases:
* Section 6.1. Gradient-dependent diffusivity case (3.1): mean curvature flow of Example 3.2, a singular reaction-source term, inhomogeneous Dirichlet boundary conditions, and known exact solution.

* Section 6.2. Gradient-independent diffusivity case (3.3): one time step of a backward Euler time discretization of the Richards equation with nonlinear degenerate advection-reaction-diffusion terms, anisotropic permeability tensor, and mixed boundary conditions.

## 6.1 Gradient-dependent diffusivity

<figure><img src="image_2.png" /><figcaption>Figure 1: [Section 6.1] The computational domain and the exact solution for the gradient-dependent diffusivity test case.</figcaption></figure>

The computational domain $\Omega$ is the square $(-1, 1)^2$ with a triangular section removed, see Figure 1. In $\Omega$, we consider, for (3.1),

$$ \sigma(\mathbf{x}, \nabla \xi) := A(|\nabla \xi|)\nabla \xi \quad \text{with} \quad A(\varrho) := 2\tau + 1/\sqrt{1 + \varrho^2}, \qquad (6.3a) $$

following Example 3.2, where $\tau$ is a parameter. Thus $\sigma_m = 2\tau$, $\sigma_M - \sigma_m = 1$, and Assumption (A1) is satisfied. We define

$$ f(\mathbf{x}, \xi) := \nu \xi - g(\mathbf{x}), \qquad (6.3b) $$

where $\nu = 10^{-2}$ is a fixed parameter. The function $g$ follows from our choice of the exact solution: Let $\mathbf{x} \in \Omega$ be represented in polar coordinates by the radius-azimuth pair $(r, \theta)$ and let $\lambda := 4/7$ be a fixed parameter. Then,

$$ u_{\text{exact}}(\mathbf{x}) := r^{\lambda} \cos(\lambda\theta), \quad \text{yielding} \quad g(\mathbf{x}) := \left[ r^{-1} \frac{(1-\lambda)(\lambda r^{\lambda-1})^3}{(1+(\lambda r^{\lambda-1})^2)^{\frac{3}{2}}} + \nu r^{\lambda} \right] \cos(\lambda\theta). \qquad (6.4) $$

Observe that $g \in H^{-1}(\Omega)$ but $g \notin L^2(\Omega)$ due to the singularity at $r = 0$. Thus, Assumption (A2) is not satisfied, but our numerical results seem to carry over even to this case. We set the Dirichlet condition $u = u_{\text{exact}}$ on $\partial\Omega$.

We consider two linearization schemes from Table 1, the Kačanov and the Zarantonello ones. In the Zarantonello case, we choose $\Lambda = 1$; this does not satisfy the theoretical requirement $\Lambda \ge \sigma_M^2/\sigma_m$ from the a priori analysis [66] but numerically leads to a fast convergence (no requirement on $\Lambda$ appears in our a posteriori theory). In fact, for the theoretical choice of $\Lambda = \sigma_M^2/\sigma_m$, the Zarantonello scheme barely converges for $\tau < 10^{-1}$. The initial guess is taken to be $u_\ell^0 = \frac{11}{10}u_{\text{exact}}$.

For $\tau = 0.01$, Figure 2 (left) plots the global effectivity indices for these schemes against the linearization iteration index $i$. The indices are between 1–2 in all cases. Figure 2 (center) shows the distribution of the local errors. As expected, the error is concentrated around the
<figure><img src="image_3.png" /><figcaption>Figure 2: [Section 6.1, $\tau = 0.01$] Comparison between errors and estimators. The Kačanov linearization (top) and the Zarantonello linearization (bottom). Global effectivity indices (6.2a) as functions of the refinement level $\ell$ and the linearization iteration $i$ (left). Elementwise errors $\|u_{\ell}^i - u_{\text{ref}}^{i+1}\|_{1, u_{\ell}^i, K}$ (center) and estimators $\eta_K^i$ (right) for $\ell = 2$ and at the linearization iteration $i = \bar{i}$ satisfying the stopping criterion (5.16) (both in logarithmic scale).</figcaption></figure>

singular re-entrant corner at the origin. It decreases radially, and its minima are situated at the other corners. The distribution of the estimators $\eta_K^i$, shown in Figure 2 (right), reproduces the same pattern. Figure 3 shows the corresponding local effectivity indices (6.2b). In almost every mesh element, these indices are between 0.75–1.5, and in all the cases they are in the order of unity. The effectivity indices around the singularity at the origin are about 1.5–2.0 in all the simulations, which we find to be particularly impressive. Note that an effectivity index less than 1 does not contradict the local efficiency estimate of Theorem 5.7 since, instead of the neighborhood $\tilde{K}$, the estimates are computed over the individual mesh elements $K \in \mathcal{T}_\ell$.

The robustness of the estimates is next exhibited in Figure 4. For a fixed mesh, it plots the global effectivity indices (6.2a) of the linearization schemes against $1/\tau$, where we consider $\tau \in [10^{-3}, 1]$. Though the Lipschitz continuity/monotonicity ratio $\sigma_M/\sigma_m = (1+2\tau)/2\tau$, expressing the strength of the nonlinearity, varies over three orders of magnitude, the effectivity indices stay essentially intact. This illustrates our robustness result from Theorem 5.7. More precisely, for the Zarantonello linearization, we know from Remark 5.9 that the variability constant $\vartheta_{\Omega,\ell,i}$ equals one, so that our estimates are robust. For the Kačanov scheme, as explained in Remark 5.12, we can first compute the variability constants $\vartheta_{\Omega,\ell,i}$ to a posteriori assess the robustness also in this case, see Figure 4 (left) for the values of $\vartheta_{\Omega,\ell,i}$. We recall that $\vartheta_{\Omega,\ell,i}$ is from (5.15) computed easily and locally, and without the knowledge of $u_{\text{ref}}^{i+1}$ or $u_{\text{exact}}$ needed to compute the (global) effectivity indices (6.2a).

Figure 5 further illustrates how the errors, the linearization estimator $\eta_{\text{lin},\Omega}^i$, and the total estimator $\eta_{\Omega}^i$ vary with the linearization iterations. Linearization error dominates the total error at the start of the iterations but becomes negligible as the iterations converge, where the discretization component becomes dominant. The appropriateness of (5.16) is evident from the plots, since the adaptive criterion is met only after 3 iterations for the Kačanov scheme, and after 7–9 iterations for the Zarantonello scheme, and there is no significant change in the total
<figure><img src="image_4.png" /><figcaption>Figure 3: [Section 6.1, $\tau = 0.01$] Local effectivity indices (6.2b) for the Kačanov (top) and the Zarantonello (bottom) linearizations at different discretization levels $\ell$ and for $i = \bar{i}$ from (5.16).</figcaption></figure>

errors beyond this point. Furthermore, from the simulations, the Kačanov scheme is found to be much faster than the Zarantonello scheme. The Zarantonello scheme, though, leads to slightly better effectivity indices, probably due to the local variability constants $\vartheta_{a,\ell,i}$ from (5.15) being 1 for this scheme (see Remark 5.9).

Finally, Figure 6 shows the relative magnitude of the components of the discretization estimator $\eta_{\text{disc},\Omega}^i$ defined in (5.13a). The flux estimator $\eta_{F,\Omega}^i$ is the largest contributor to $\eta_{\text{disc},\Omega}^i$ followed by the data-oscillation estimator $\eta_{\text{osc},\Omega}^i$. Since, in Table 1, $L^i = \nu$ for the Kačanov scheme and $L^i = 0$ for the Zarantonello scheme, we have that $\eta_{\text{quad},S,\Omega}^i = 0$ in both cases, whereas $\eta_{\text{quad},F,\Omega}^i < 10^{-5}$. With $L^i = 0$, the Zarantonello scheme is a pure diffusion linearization, see Remarks 4.2 and 5.5, so that $\eta_{S,\Omega}^i = 0$.

## 6.2 Gradient-independent diffusivity: the Richards equation

We choose here $\Omega$ as the unit square, $\Omega = (0,1) \times (0,1)$. The backward Euler time-discretized Richards equation is characterized here by the following choices of the functions in (3.3):

$$ 
\begin{gathered}
f(\boldsymbol{x}, \boldsymbol{\xi}) = S(\boldsymbol{\xi}) - S(\bar{u}(\boldsymbol{x})), \quad \mathcal{D}(\boldsymbol{x}, \boldsymbol{\xi}) = \kappa(S(\boldsymbol{\xi})), \quad \boldsymbol{q}(\boldsymbol{x}, \boldsymbol{\xi}) = -\kappa(S(\boldsymbol{\xi})) \boldsymbol{g}, \\
\bar{K} = \begin{bmatrix} 1 & 0.2 \\ 0.2 & 1 \end{bmatrix}, \quad \boldsymbol{g} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \quad \tau \in [10^{-3}, 1].
\end{gathered}
\tag{6.5} $$

Here $\boldsymbol{g}$ is the gravity vector, and the functions $S: \mathbb{R} \to (0,1]$ and $\kappa: [0,1] \to [0,1]$, known as the saturation and relative permeability functions, respectively, are modeled using the van Genuchten [62] parametrization, expressed in a nondimensional setting for $\lambda \in (0,1)$ as

$$ S(\xi) := \left(1 + (2 - \xi)^{\frac{1}{1-\lambda}}\right)^{-\lambda}, \quad \kappa(s) := \sqrt{s} \left(1 - (1 - s^{\frac{1}{\lambda}})^{\lambda}\right)^2. \tag{6.6} $$
<figure><img src="image_5.png" /><figcaption>Figure 4: [Section 6.1, $\tau$ varies] Global effectivity indices (6.2a) vs. $1/\tau$ for the Kačanov (left) and Zarantonello (right) linearizations at different discretization levels $\ell$ for $i = \bar{i}$ from (5.16). For the Kačanov scheme, the variability constant $\vartheta_{\Omega,\ell,i}$ from (5.15) is also plotted for $\ell = 2$. For Zarantonello linearization, this constant is always 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: [Section 6.1, $\tau = 0.01$] Total error $\|\mathcal{R}(u_{\ell}^i)\|_{-1, u_{\ell}^i}$ approximated using (6.1), linearization error $\eta_{\text{lin},\Omega}^i$, and the total estimator $\eta_{\Omega}^i$ against iteration index $i$ in a log$_{10}$-scale: the Kačanov scheme (top), the Zarantonello scheme (bottom). The vertical black dashed line on each plot indicates the stopping criterion $i = \bar{i}$ from (5.16).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: [Section 6.1, $\tau = 0.01$] The components (5.12) of the discretization estimator $\eta_{\text{disc},\Omega}^i$ vs. iteration $i$ for $\tau = 0.01$. For both schemes, $\eta_{\text{quad,S},\Omega}^i = 0$, $\eta_{\text{quad,F},\Omega}^i < 10^{-5}$, and for the Zarantonello scheme, $\eta_{S,\Omega}^i = 0$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: [Section 6.2] The domain $\Omega = (0,1)^2$ with the boundary conditions and reference solutions of (1.1) for $\tau = 1$ (left) and $\tau = 0.01$ (right).</figcaption></figure>

More details on the Richards equation can be found in [53] and the references therein. For the rest of this paper, we will use $\lambda = 0.5$. The mixed (Dirichlet and no-flux) boundary condition used for the problem is shown in Figure 7 with $u_{\text{in}} = 1$ and $u_{\text{out}} = -3$. Finally, in the context of the time-dependent Richards equation, the quantity $\bar{u}$ represents the solution of the previous time step. Here, it is arbitrarily taken to be the solution of the Darcy problem $\nabla \cdot [\bar{K}(\nabla \bar{u} - g)] = 0$ in $\Omega$, completed with the mixed boundary condition mentioned above. It is straightforward to verify that the Richards equation satisfies the conditions stated in (B1)–(B3). Furthermore, it satisfies (B4) since in order for the constant $\gamma > 0$ to exist, one requires by Cauchy’s mean value theorem that $(\kappa'(S(\xi))^2 S'(\xi) / \kappa(S(\xi))$ is bounded for all $\xi \in \mathbb{R}$, a requirement that is satisfied for $\lambda = 0.5$ ($\gamma \approx 1.2$ in this case).

<figure><img src="image_9.png" /><figcaption>Figure 8: [Section 6.2, $\tau = 1$ and $\tau = 0.01$] Global effectivity indices (6.2) for different linearization schemes, $\tau$, $\ell$, and $i$ values. Results for $\tau = 1$ (top) and $\tau = 0.01$ (bottom). The Picard scheme (left), the *M*-scheme (center), and the *L*-scheme (right).</figcaption></figure>

We take the initial guess for the iterations as $u_{\ell}^0 = 0$ and show results for numerical schemes mentioned in Table 2. We fix the parameters for *L*- and *M*-schemes as $L = 1/6$ and $M = 0.1$ respectively. The *M*-scheme is chosen to depict the base numerical results since it is a hybrid
between the L-scheme and the Picard (fixed point) iterations [52]. Global effectivity indices, computed using (6.2a), are plotted in Figure 8 for the different schemes, at all discretization levels, different time-step lengths $\tau$, and iterations. It is seen that the effectivity indices stay above 1 in all cases, thus validating the guaranteed upper bound (5.14a). As the linearization iterations converge, these indices reach a stable value as expected.

To inspect the robustness of our estimators with respect to the strength of the nonlinearities, we in Figure 9 plot the effectivity indices against the Lipschitz/monotonicity constant of the operator $\mathcal{R}$, represented in this case by $1/\tau$. Despite the parameter $\tau$ varying between 1 and $10^{-3}$, thus approaching the singularly perturbed regime, the effectivity indices stay bounded between 1.5 and 2.6. From Table 2 and (5.15), the variability constant $\vartheta_{\Omega,\ell,i}$ essentially depends on the patchwise variations of $\mathcal{D}(x, u_{\ell}^i) = \kappa(S(u_{\ell}^i))$. Scaled by half, we plot $\vartheta_{\Omega,\ell,i}$ for all the schemes on the second mesh level $\ell = 2$. We observe uniformly moderate values, which immediately ensure the observed robustness, see Remark 5.12 (recall that $\vartheta_{\Omega,\ell,i}$ is easily and locally computable, whereas we can only compute the (global) effectivity indices (6.2a) knowing/approximating the exact solution).

<figure><img src="image_10.png" /><figcaption>Figure 9: [Section 6.2, $\tau$ varies] Global effectivity indices (6.2a) vs. $1/\tau$ for the Picard (left), M- (center), and L-schemes (right) at different discretization levels $\ell$ for $i = \bar{i}$ from (5.16). The variability constant $\vartheta_{\Omega,\ell,i}$ from (5.15) is also plotted (scaled by half and for $\ell = 2$).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: [Section 6.2, $\tau = 1$ and $\tau = 0.01$] The elementwise error approximated by $\|u_{\ell}^i - u_{\text{ref}}^{i+1}\|_{1,u_{\ell}^i,K}$ vs. the elementwise estimators $\eta_K^i$ for $\ell = 2$ for the M-scheme at $i = \bar{i}$ from (5.16); $\tau = 1$ (left) and $\tau = 0.01$ (right) (all in logarithmic scale).</figcaption></figure>

Figure 10 presents the distribution of the error for the M-scheme at $\ell = 2$ compared to the distribution of $\eta_K^i$. Using the maximum principle, we obtain from (6.6) that $\min_{\Omega} \kappa(S(u)) = \kappa(S(u_{\text{out}})) \approx 5 \times 10^{-4}$, whereas $\max_{\Omega} \kappa(S(u)) = \kappa(S(u_{\text{in}})) = 1$. As a result, the elementwise errors and estimators are seen to vary 2–4 orders of magnitude. Nevertheless, the plots suggest that the estimator $\eta_K^i$ captures both the magnitude and the distribution of the error accurately. This is quantified by the local effectivity indices (6.2b) plotted in Figure 11. In most of the domain, the effectivity index is between 1 and 2.5, and in all cases it is in the order of unity. For the other linearization schemes, the local effectivity indices also stay within this range, despite the estimator $\eta_{\Omega}^i$ varying by about 2 to 4 orders of magnitude, $\tau$ varying 3 orders of magnitude, and the diffusion coefficient $\kappa \circ S$ in (6.5)–(6.6) varying by 3 to 4 orders of magnitude. Hence, we consider the match to be excellent.
<figure><img src="image_12.png" /><figcaption>Figure 11: [Section 6.2, $\tau = 1$ and $\tau = 0.01$] Local effectivity indices (6.2b) for the M-scheme at discretization levels $\ell = 1, 2, 4$ for $i = \bar{i}$ from (5.16); $\tau = 1$ (top) and $\tau = 0.01$ (bottom).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: [Section 6.2, $\tau = 1$ and $\tau = 0.01$] Total error $|||\mathcal{R}(u_{\ell}^i)|||_{-1, u_{\ell}^i}$ approximated using (6.1), linearization error $\eta_{\text{lin},\Omega}^i$, and total estimator $\eta_{\Omega}^i$ against iteration index $i$ in a log$_{10}$-scale. Fixed $\ell = 2$. The vertical black dashed line indicates the stopping criterion $i = \bar{i}$ from (5.16). Results are for $\tau = 1$ (top) and $\tau = 0.01$ (bottom). The Picard scheme (left), the M-scheme (center), and the L-scheme (right).</figcaption></figure>

Next, we investigate how the error, the linearization estimator $\eta_{\text{lin},\Omega}^i$, and the total estimator $\eta_{\Omega}^i$ vary with the iteration $i$ for different linearization schemes. This is shown in Figure 12. At the beginning, the dominant error component is the linearization error, which is roughly equal to the total error. However, as the iterations advance, the total error is first dominated by the
discretization error, and then dictated by it. The results also clearly portray the known trends of the linearization schemes [50, 52], i.e., the L-scheme converges more slowly as $\tau$ is decreased, whereas the Picard and the M-schemes speed up, and the M-scheme is faster than the Picard scheme particularly for larger values of $\tau$. Figure 12 further illustrates the efficiency of (5.16), since the iterations are stopped after they cease to have any significant effect on the quality of the solutions.

<figure><img src="image_14.png" /><figcaption>Figure 13: [Section 6.2, $\tau = 1$ and $\tau = 0.01$] The relative magnitudes of components of the discretization estimator $\eta_{\text{disc},\Omega}^i$ from Section 5.2 plotted against iteration $i$ for the M-scheme at discretization level $\ell = 2$; $\tau = 1$ (left) and $\tau = 0.01$ (right).</figcaption></figure>

Finally, Figure 13 shows the relative magnitude of the components of the discretization estimator $\eta_{\text{disc},\Omega}^i$, see (5.13a). The flux estimator $\eta_{F,\Omega}^i$ is the dominant component in all cases, followed in order by the potential estimator $\eta_{S,\Omega}^i$, the data-oscillation estimator $\eta_{\text{osc},\Omega}^i$, and the quadrature estimator $\eta_{\text{quad},F,\Omega}$ (we also plot $\eta_{\text{quad},S,\Omega}$).

# 7 Well-posedness of the models (3.1) and (3.3)

In this section, we collect the existence and uniqueness proofs for the models (3.1) and (3.3).

## 7.1 Well-posedness of the gradient-dependent diffusivity model (3.1)

**Proposition 7.1** (Monotonicity and Lipschitz continuity of $\mathcal{R}$ when defined by (3.1), existence and uniqueness of the weak solution of (1.1)). Assume the expression (3.1) for $\mathcal{R}: H_0^1(\Omega) \to H^{-1}(\Omega)$ and (A1)-(A2). Then, for all $\xi, \zeta \in H_0^1(\Omega)$,

$$ \sigma_m \|\nabla(\xi - \zeta)\| \le \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \varphi \rangle}{\|\nabla \varphi\|} \le (\sigma_M + (C_\Omega h_\Omega)^2 f_M) \|\nabla(\xi - \zeta)\|, \quad (7.1) $$

and there exists a unique solution $u \in H_0^1(\Omega)$ of (1.1).
Proof. The operator $\mathcal{R}$ from (3.1) is monotone and continuous since

$$
\begin{align*}
\langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \xi - \zeta \rangle &\stackrel{(3.1)}{=} (\sigma(\cdot, \nabla\xi) - \sigma(\cdot, \nabla\zeta), \nabla(\xi - \zeta)) + \langle f(\cdot, \xi) - f(\cdot, \zeta), \xi - \zeta \rangle \\
&\stackrel{(A2)}{\ge} (\sigma(\cdot, \nabla\xi) - \sigma(\cdot, \nabla\zeta), \nabla(\xi - \zeta)) \stackrel{(A1)}{\ge} \sigma_m \|\nabla(\xi - \zeta)\|^2, \\
\sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \varphi \rangle}{\|\nabla \varphi\|} &\stackrel{(A1)}{\le} \sigma_M \|\nabla(\xi - \zeta)\| + \sup_{\varphi \in H_0^1(\Omega)} \frac{\|f(\cdot, \xi) - f(\cdot, \zeta)\|\|\varphi\|}{\|\nabla \varphi\|} \\
&\stackrel{(2.1)}{\le} \sigma_M \|\nabla(\xi - \zeta)\| + C_\Omega h_\Omega \|f(\cdot, \xi) - f(\cdot, \zeta)\| \\
&\stackrel{(A2)}{\le} \sigma_M \|\nabla(\xi - \zeta)\| + C_\Omega h_\Omega f_M \|\xi - \zeta\| \\
&\stackrel{(2.1)}{\le} (\sigma_M + (C_\Omega h_\Omega)^2 f_M) \|\nabla(\xi - \zeta)\|.
\end{align*}
$$

The existence of $u \in H_0^1(\Omega)$ solving (1.1) then follows from the nonlinear Lax-Milgram theorem [67, Chapter 2] due to the monotonicity and Lipschitz continuity of $\mathcal{R}$. $\square$

## 7.2 Well-posedness of the gradient-independent diffusivity model (3.3)

**Proposition 7.2** (Monotonicity and Lipschitz continuity of $\mathcal{R}$ when defined by (3.3), existence and uniqueness of the weak solution of (1.1)). Assume the expression (3.3) for $\mathcal{R}: H_0^1(\Omega) \to H^{-1}(\Omega)$ and (B1)-(B5) with $0 < \tau < (4/\gamma)$. Moreover, let the Kirchhoff transform function $\mathcal{K}: \Omega \times \mathbb{R} \to \mathbb{R}$, defined as

$$ \mathcal{K}(\mathbf{x}, \xi) := \int_{0}^{\xi} \mathcal{D}(\mathbf{x}, \varrho) \, d\varrho, \qquad (7.2) $$

satisfy $\lim_{\xi \to \pm\infty} \mathcal{K}(\mathbf{x}, \xi) = \pm\infty$ for all $\mathbf{x} \in \Omega$. Then, for all $\xi, \zeta \in H_0^1(\Omega)$,

$$ \|\nabla(\mathcal{K}(\cdot, \xi) - \mathcal{K}(\cdot, \zeta))\| \lesssim \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \varphi \rangle}{\|\nabla \varphi\|} \lesssim \|\nabla(\mathcal{K}(\cdot, \xi) - \mathcal{K}(\cdot, \zeta))\|, \quad (7.3) $$

where the hidden constants depend only on $\Omega$ and the constants $\tau, \gamma, \mathcal{D}_M, f_M, K_M, K_M$, and $q_M$. Moreover, there exists a unique solution $u \in L^\infty(\Omega) \cap H_0^1(\Omega)$ of (1.1).

**Proof.** To prove monotonicity of $\mathcal{R}$, we introduce $\mathbf{v}(\mathbf{x}, \xi) := \mathbf{q}(\mathbf{x}, \xi) - \sum_{j=1}^d e_j \int_0^\xi \partial_{x_j} \mathcal{D}(\mathbf{x}, \varrho) \, d\varrho$ to rewrite (3.3) using (7.2) and $\nabla \mathcal{K}(\mathbf{x}, \xi) = \mathcal{D}(\mathbf{x}, \xi) \nabla \xi + \sum_{j=1}^d e_j \int_0^\xi \partial_{x_j} \mathcal{D}(\mathbf{x}, \varrho) \, d\varrho$ as

$$ \langle \mathcal{R}(\xi), \varphi \rangle = \tau(\bar{\mathbf{K}}(\cdot)(\nabla \mathcal{K}(\cdot, \xi) + \mathbf{v}(\cdot, \xi)), \nabla \varphi) + \langle f(\cdot, \xi), \varphi \rangle. \quad (7.4) $$

Dropping the $\mathbf{x}$-dependence from the functions $f, \bar{\mathbf{K}}, \mathcal{K}$, and $\mathbf{v}$ for simplicity, we have for any $\xi, \zeta \in H_0^1(\Omega)$,

$$ \langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \varphi \rangle = \langle f(\xi) - f(\zeta), \varphi \rangle + \tau(\bar{\mathbf{K}}(\nabla(\mathcal{K}(\xi) - \mathcal{K}(\zeta)) + \mathbf{v}(\xi) - \mathbf{v}(\zeta)), \nabla \varphi). \quad (7.5) $$

For obtaining the monotonicity bound, we test (7.5) with $\varphi = \mathcal{K}(\xi) - \mathcal{K}(\zeta) = \int_\zeta^\xi \mathcal{D} \stackrel{(B1)}{\in} H_0^1(\Omega)$ and use the Young inequality (2.2) with $\rho = 2/(\tau\gamma)$, where $\gamma > 0$ is defined in (B4):

$$
\begin{align}
\langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \mathcal{K}(\xi) - \mathcal{K}(\zeta) \rangle \nonumber \\
&\stackrel{(7.5)}{=} \langle f(\xi) - f(\zeta), \int_{\zeta}^{\xi} \mathcal{D} \rangle + \tau \|\bar{\mathbf{K}}^{\frac{1}{2}} \nabla(\mathcal{K}(\xi) - \mathcal{K}(\zeta))\|^2 + \tau(\bar{\mathbf{K}}(\mathbf{v}(\xi) - \mathbf{v}(\zeta)), \nabla(\mathcal{K}(\xi) - \mathcal{K}(\zeta))) \nonumber \\
&\stackrel{(2.2)}{\ge} \langle f(\xi) - f(\zeta), \int_{\zeta}^{\xi} \mathcal{D} \rangle + \left(1 - \frac{\tau\gamma}{4}\right) \tau \|\bar{\mathbf{K}}^{\frac{1}{2}} \nabla(\mathcal{K}(\xi) - \mathcal{K}(\zeta))\|^2 - \gamma^{-1} \|\bar{\mathbf{K}}^{\frac{1}{2}}(\mathbf{v}(\xi) - \mathbf{v}(\zeta))\|^2 \tag{7.6} \\
&\stackrel{(B4)}{\ge} \left(1 - \frac{\tau\gamma}{4}\right) \tau \|\bar{\mathbf{K}}^{\frac{1}{2}} \nabla(\mathcal{K}(\xi) - \mathcal{K}(\zeta))\|^2 \stackrel{(B3)}{\ge} \left(1 - \frac{\tau\gamma}{4}\right) \tau K_M \|\nabla(\mathcal{K}(\xi) - \mathcal{K}(\zeta))\|^2. \nonumber
\end{align}
$$
Hence, if $\tau < 4/\gamma$, then we have the lower bound by using $\varphi = K(\xi) - K(\zeta)$ in (7.3).

For obtaining the Lipschitz-continuity result, we first define for some $\varepsilon \in (0, D_M)$ (with $D_M > 0$ from (B1)) and $\xi \in \mathbb{R}$,

$$ D_{\varepsilon}(\xi) := \max(\varepsilon, D(\xi)) \quad \text{and} \quad K_{\varepsilon}(\xi) := \int_{0}^{\xi} D_{\varepsilon}(\varrho) d\varrho. \qquad (7.7a) $$

Since, for $\zeta \in \mathbb{R}$,

$$ |K(\xi) - K(\zeta)| = \left| \int_{\zeta}^{\xi} D \right| \le \left| \int_{\zeta}^{\xi} D_{\varepsilon} \right| = |K_{\varepsilon}(\xi) - K_{\varepsilon}(\zeta)| \le \|D_{\varepsilon}\|_{L^{\infty}} |\xi - \zeta|, \qquad (7.7b) $$

we have the following properties of the regularized function $K_\varepsilon$:

$$ \max(\varepsilon|\xi - \zeta|, |K(\xi) - K(\zeta)|) \le |K_{\varepsilon}(\xi) - K_{\varepsilon}(\zeta)| \stackrel{(B1)}{\le} D_M |\xi - \zeta|. \qquad (7.7c) $$

Then, substituting $K$ by $K_\varepsilon$ in (7.5), the operator $\mathcal{R}$ still remains monotone following the steps of (7.6). Moreover, repeatedly using the Poincaré inequality (2.1), we obtain

$$
\begin{align*}
& \sup_{\substack{\varphi \in H_0^1(\Omega), \\ \|\nabla \varphi\| = 1}} \langle \mathcal{R}(\xi) - \mathcal{R}(\zeta), \varphi \rangle \\
& \stackrel{(B3)}{\le} \sup_{\substack{\varphi \in H_0^1(\Omega), \\ \|\nabla \varphi\| = 1}} \|f(\xi) - f(\zeta)\| \|\varphi\| + \tau K_M \|\nabla(\mathcal{K}_\varepsilon(\xi) - \mathcal{K}_\varepsilon(\zeta))\| \\
& \phantom{\stackrel{(B3)}{\le}} + \tau K_M^{\frac{1}{2}} \|\bar{\mathbf{K}}^{\frac{1}{2}}(\mathbf{v}(\xi) - \mathbf{v}(\zeta))\| \\
& \stackrel{(2.1),(B4)}{\le} C_\Omega h_\Omega \|f(\xi) - f(\zeta)\| + \tau K_M \|\nabla(\mathcal{K}_\varepsilon(\xi) - \mathcal{K}_\varepsilon(\zeta))\| \\
& \phantom{\stackrel{(2.1),(B4)}{\le}} + \tau(\gamma K_M)^{\frac{1}{2}} \|f(\xi) - f(\zeta)\|^{\frac{1}{2}} \|\mathcal{K}(\xi) - \mathcal{K}(\zeta)\|^{\frac{1}{2}} \\
& \stackrel{(B1),(B2),(7.7b),(2.2)}{\lesssim} \|\xi - \zeta\| + \tau K_M \|\nabla(\mathcal{K}_\varepsilon(\xi) - \mathcal{K}_\varepsilon(\zeta))\| \\
& \stackrel{(7.7c)}{\le} \frac{1}{\varepsilon} \|\mathcal{K}_\varepsilon(\xi) - \mathcal{K}_\varepsilon(\zeta)\| + \tau K_M \|\nabla(\mathcal{K}_\varepsilon(\xi) - \mathcal{K}_\varepsilon(\zeta))\| \\
& \stackrel{(2.1)}{\le} \left( \frac{C_\Omega h_\Omega}{\varepsilon} + \tau K_M \right) \|\nabla(\mathcal{K}_\varepsilon(\xi) - \mathcal{K}_\varepsilon(\zeta))\|.
\end{align*}
$$

Then, the existence of $u_\varepsilon$ with $K_\varepsilon(u_\varepsilon) \in H_0^1(\Omega)$ solving (1.1) follows from the monotonicity and continuity of $\mathcal{R}$ using the nonlinear Lax-Milgram theorem [67, Chapter 2] (see also [68, Chapter 25]). Moreover, $K_\varepsilon(u_\varepsilon) \in L^\infty(\Omega)$, which is evident from Theorem 13.1 of [47, Chapter 3], since $\mathbf{q}$, $\partial x_j D$, $f$ are all bounded in (7.4). In particular, since the $L^\infty$ bounds on these functions do not depend on $\varepsilon$, the bound on $\|K_\varepsilon(u_\varepsilon)\|_{L^\infty(\Omega)}$ is also independent of $\varepsilon$. Then, from $|K_\varepsilon(\xi)| \ge |K(\xi)|$ (as evident from putting $\zeta = 0$ in (7.7c)), we have that $\|K(u_\varepsilon)\|_{L^\infty(\Omega)}$ is bounded irrespective of $\varepsilon$. The assumption $K(\xi) \to \pm\infty$ as $\xi \to \pm\infty$ further implies that there exists a constant $R \in (0, \infty)$ independent of $\varepsilon > 0$ such that $\|u_\varepsilon\|_{L^\infty(\Omega)} \le R$. Define $\min_{|v| \le R} D(v) := \varepsilon_0 > 0$, the positivity following due to (B1) since $\{|v| \le R\}$ is compact. Then, for any $0 < \varepsilon < \varepsilon_0$, we have that $u_\varepsilon$ also solves the unregularized problem as $D(u_\varepsilon) \ge \varepsilon_0 > \varepsilon$, and thus, $K(u_\varepsilon) = K_\varepsilon(u_\varepsilon)$ from (7.7a). Moreover, $u = u_\varepsilon \in H_0^1(\Omega)$ (as $(K_\varepsilon)^{-1}$ is Lipschitz from (7.7c)) and it is the unique solution due to the monotonicity of the operator $\mathcal{R}$, see (7.6). $\square$

# 8 A posteriori error estimates of Theorem 5.7

We show in this section how Theorem 5.3 together with the results on a posteriori error analysis of linear reaction-diffusion problems lead to the a posteriori error estimates of the nonlinear problem (1.1) given in Theorem 5.7.
## 8.1 Equilibrated flux and potential of Definition 5.4, proof of Lemma 5.6

Let us first give all the details for the construction of Definition 5.4.

### 8.1.1 Orthogonal projections

Let the standard $L^2$-orthogonal projections $\Pi_{\ell,p} : L^2(\Omega) \to \mathcal{P}_p(\mathcal{T}_\ell)$ and $\Pi_{\ell,p-1}^{\text{RT}} : L^2(\Omega; \mathbb{R}^d) \to \mathcal{RT}_{p-1}(\mathcal{T}_\ell; \mathbb{R}^d)$ be defined as

$$ \text{for } \zeta \in L^2(\Omega), \quad (\Pi_{\ell,p}\zeta, v_\ell) = (\zeta, v_\ell) \quad \forall v_\ell \in \mathcal{P}_p(\mathcal{T}_\ell), \quad (8.1a) $$

$$ \text{for } \zeta \in L^2(\Omega; \mathbb{R}^d), \quad (\Pi_{\ell,p-1}^{\text{RT}}\zeta, v_\ell) = (\zeta, v_\ell) \quad \forall v_\ell \in \mathcal{RT}_{p-1}(\mathcal{T}_\ell; \mathbb{R}^d). \quad (8.1b) $$

Note that (8.1) can be equivalently rewritten elementwise. We will also need the $L^i$-weighted projection $\Pi_{\ell,p,L^i} : \mathcal{P}_{p+1}(\mathcal{T}_\ell) \to \mathcal{P}_p(\mathcal{T}_\ell)$

$$ \text{for } \zeta_\ell \in \mathcal{P}_{p+1}(\mathcal{T}_\ell), \quad (L^i \Pi_{\ell,p,L^i} \zeta_\ell, v_\ell) = (L^i \zeta_\ell, v_\ell) \quad \forall v_\ell \in \mathcal{P}_p(\mathcal{T}_\ell), \quad (8.1c) $$

defined again elementwise and meaning, for all $K \in \mathcal{T}_\ell$,

$$ \Pi_{\ell,p,L^i} \zeta_\ell|_K = \zeta_\ell|_K \quad \text{when } L^i \text{ is zero on } K, \\ (L^i \Pi_{\ell,p,L^i} \zeta_\ell, v_\ell)_K = (L^i \zeta_\ell, v_\ell)_K \quad \forall v_\ell \in \mathcal{P}_p(K) \text{ otherwise.} $$

This last projection only serves to describe data oscillations. When $L^i$ is elementwise constant, then $\Pi_{\ell,p,L^i} = \Pi_{\ell,p}$.

### 8.1.2 Trace and inverse inequalities

For all $K \in \mathcal{T}_\ell$ there exists a constant $C_{\text{Tr}} > 0$ depending only on the space dimension $d \ge 1$ and shape-regularity of the mesh $\mathcal{T}_\ell$ such that the there holds

$$ \text{Trace inequality: } \|v\|_{\partial K} \le C_{\text{Tr}} \|\nabla v\|_{K}^{\frac{1}{2}} \|v\|_{K}^{\frac{1}{2}} \quad \forall v \in H^1(K) \text{ with } \int_K v = 0. \quad (8.2a) $$

Similarly, for any $v \in \mathcal{RT}_p(K)$, there exist constants $C_{\text{inv},p}$, $C_{\text{inv},p,\partial} > 0$ which depend only on $d$, the shape-regularity of the mesh, and the polynomial degree $p$ such that

$$ \text{Inverse inequalities: } \begin{cases} h_K \|\nabla \cdot v\|_K \le C_{\text{inv},p} \|v\|_K, \\ h_K^{\frac{1}{2}} \|v \cdot n\|_{\partial K} \le C_{\text{inv},p,\partial} \|v\|_K. \end{cases} \quad (8.2b) $$

Possible expressions of the constants in (8.2) are provided in Section 2.2 of [60], along with the proofs of these statements.

### 8.1.3 Weights

As shown in [60, 64], in order to obtain robust reliability and efficiency estimates for singularly perturbed linear reaction-diffusion problems (for large values of $L^i$), weights need to be introduced. Recalling (4.5), we define a piecewise constant function $w^i$

$$ w^i|_K := \min \left( 1, \frac{C_*}{h_K^{\frac{1}{2}}} \left( \frac{a_m^i|_K}{L_m^i|_K} \right)^{\frac{1}{4}} \right), \quad K \in \mathcal{T}_\ell, \quad (8.3a) $$

where, with reference to constants introduced in (8.2),

$$ C_* := \frac{1}{\sqrt{2}} \left( \frac{1}{\sqrt{\pi}} C_{\text{inv},p} + C_{\text{Tr}} C_{\text{inv},p,\partial} \right). $$

To treat data oscillation terms, the usual cut-off factor [64] for $K \in \mathcal{T}_\ell$ is

$$ \tilde{w}_K^i := \min \left( \frac{1}{(L_m^i|_K)^{\frac{1}{2}}}, \frac{h_K}{\pi (a_m^i|_K)^{\frac{1}{2}}} \right). \quad (8.3b) $$
### 8.1.4 Shorhand notations

It will be useful to use below the shortened notation

$$ \tau^{i+1} := \mathbf{a}^i \nabla u_{\ell}^{i+1} + \mathcal{F}^i. \qquad (8.4) $$

We will also employ, for $K \in \mathcal{T}_\ell$, as in (4.19a),

$$ \|\varphi\|_{i,K} := ((L^i \varphi, \varphi)_K + (\mathbf{a}^i \nabla \varphi, \nabla \varphi)_K)^{\frac{1}{2}}. \qquad (8.5) $$

Finally, for a mesh element $K \in \mathcal{T}_\ell$, we denote its set of vertices by $\mathcal{V}_K \subset \mathcal{V}_\ell$.

### 8.1.5 Definition 5.4

The existence and uniqueness of $\sigma_a^{i+1}$ and $\phi_a^{i+1}$ from Definition 5.4 are standard following [7], see also the discussion after [60, Definition 2.3]. Problem (5.7) is of reaction-diffusion type when $L^i|_K > 0$ for at least one $K \subset \omega_a$. When $L^i|_{\omega_a} = 0$, it reduces to a pure diffusion equilibration, see Remarks 4.2 and 5.5. Note that, recalling the notation (8.4), the source term in (5.7) satisfies

$$ -(\Pi_{\ell,p}(\psi_a S^i) + \nabla \psi_a \cdot \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, 1)_{\omega_a} \stackrel{(8.1)}{=} -(S^i, \psi_a)_{\omega_a} - (\tau^{i+1}, \nabla \psi_a)_{\omega_a} \stackrel{(4.9)}{=} (L^i u_{\ell}^{i+1}, \psi_a)_{\omega_a}. \qquad (8.6) $$

Let $L^i|_{\omega_a} = 0$. Then, when $\mathbf{a} \in \mathcal{V}_\ell^{\text{int}}$ is an interior vertex, (5.6b) leads to

$$ (\nabla \cdot \sigma_a^{i+1}, 1)_{\omega_a} = (\sigma_a^{i+1} \cdot \mathbf{n}, 1)_{\partial \omega_a} = 0, \qquad (8.7) $$

and we have the necessary Neumann compatibility condition, the zero mean value of the source term in (5.7), from (8.6). Here, actually, $\phi_a^{i+1}$ becomes undefined, and we tacitly replace it by $\Pi_{\ell,p}(\psi_a u_{\ell}^{i+1})$, so that $\phi_a^{i+1} = u_{\ell}^{i+1}$. Anyhow, shall $L^i$ be zero everywhere, this is consistent since $\eta_{S,\omega}^i$ of (5.12b) features $(L^i)^{\frac{1}{2}}(u_{\ell}^{i+1} - \phi_a^{i+1}) = 0$.

Practically, $\sigma_a^{i+1}$ and $\phi_a^{i+1}$ are computed by solving the following Euler-Lagrange equations: find $(\sigma_a^{i+1}, \phi_a^{i+1}) \in \mathcal{V}_a \times \mathcal{Q}_a$ such that for all $(v_\ell, q_\ell) \in \mathcal{V}_a \times \mathcal{Q}_a$,

$$ \begin{aligned} (\sigma_a^{i+1}, (\mathbf{w}^i)^2 (a_m^i)^{-1} \mathbf{v}_\ell)_{\omega_a} - (\phi_a^{i+1}, \nabla \cdot \mathbf{v}_\ell)_{\omega_a} &= -(\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, (\mathbf{w}^i)^2 (a_m^i)^{-1} \mathbf{v}_\ell)_{\omega_a} \\ &\quad - (\psi_a u_\ell^{i+1}, \nabla \cdot \mathbf{v}_\ell)_{\omega_a}, \\ (\nabla \cdot \sigma_a^{i+1}, q_\ell)_{\omega_a} + (L^i \phi_a^{i+1}, q_\ell)_{\omega_a} &= -(\psi_a S^i + \nabla \psi_a \cdot \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, q_\ell)_{\omega_a}. \end{aligned} $$

### 8.1.6 Proof of Lemma 5.6

We have discussed the existence and uniqueness of $\sigma_a^{i+1}$ and $\phi_a^{i+1}$ in Section 8.1.5 above. The inclusions (5.9) then follow immediately from (5.8) together with the no-flow boundary conditions imposed in (5.6b).

For any vertex $\mathbf{a} \in \mathcal{V}$, the minimization constraint in (5.7) and property (8.6) yield

$$ (\nabla \cdot \sigma_a^{i+1}, 1)_{\omega_a} + (L^i \phi_a^{i+1}, 1)_{\omega_a} \stackrel{(5.7)}{=} -(\Pi_{\ell,p}(\psi_a S^i) + \nabla \psi_a \cdot \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, 1)_{\omega_a} \stackrel{(8.6)}{=} (L^i u_{\ell}^{i+1}, \psi_a)_{\omega_a}. \qquad (8.8) $$

Thus, using (8.7) for an interior vertex $\mathbf{a} \in \mathcal{V}_\ell^{\text{int}}$, we conclude the mean value property (5.10).

Finally, the minimization constraint in (5.7) together with the partition of unity

$$ \sum_{a \in \mathcal{V}_\ell} \psi_a = 1 \qquad (8.9) $$

and (5.8) imply (5.11).
## 8.2 Proof of global reliability estimate in Theorem 5.7

Recalling Theorem 5.3, we need to show that

$$ \left\| \mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1}) \right\|_{-1, u_{\ell}^{i}}^{2} \leq \sum_{K \in \mathcal{T}_{\ell}} (\eta_{\text{F},K}^{i} + \eta_{\text{S},K}^{i} + \eta_{\text{osc},K}^{i} + \eta_{\text{quad,F},K}^{i})^{2}. \quad (8.10) $$

Recall the definition of $\mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}$ from (5.1) and the norms $\|\cdot\|_{\pm 1, u_{\ell}^{i}}$ from (4.15), (4.18), and let $\varphi \in H_{0}^{1}(\Omega)$ be fixed. Using (5.2), as well as (5.9), the Green theorem, and (5.11), we have

$$ \begin{aligned} \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1}), \varphi \rangle &\stackrel{(5.2)}{=} (L^{i} u_{\ell}^{i+1} - L^{i} \phi_{\ell}^{i+1}, \varphi) + (\tau^{i+1} + \sigma_{\ell}^{i+1}, \nabla \varphi) + (S^{i} + L^{i} \phi_{\ell}^{i+1} + \nabla \cdot \sigma_{\ell}^{i+1}, \varphi) \\ &\stackrel{(5.11)}{=} (L^{i}(u_{\ell}^{i+1} - \phi_{\ell}^{i+1}), \varphi) + (\Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + \sigma_{\ell}^{i+1}, \nabla \varphi) \\ &\qquad + (\tau^{i+1} - \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, \nabla \varphi) + (S^{i} - \Pi_{\ell,p} S^{i}, \varphi) + (L^{i} \phi_{\ell}^{i+1} - \Pi_{\ell,p} (L^{i} \phi_{\ell}^{i+1}), \varphi) \\ &= \sum_{K \in \mathcal{T}_{\ell}} [(L^{i}(u_{\ell}^{i+1} - \phi_{\ell}^{i+1}), \varphi)_{K} + (\Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + \sigma_{\ell}^{i+1}, \nabla \varphi)_{K} \\ &\qquad + ((\mathbb{I} - \Pi_{\ell,p-1}^{\text{RT}}) \tau^{i+1}, \nabla \varphi)_{K} + ((I - \Pi_{\ell,p}) (S^{i} + L^{i} \phi_{\ell}^{i+1}), \varphi)_{K}]. \end{aligned} \quad (8.11) $$

We now follow the proof of [60, Theorem 3.1]. Let $K \in \mathcal{T}_{\ell}$ be fixed. For the third term of (8.11), one directly obtains, using the Cauchy–Schwarz inequality,

$$ ((\mathbb{I} - \Pi_{\ell,p-1}^{\text{RT}}) \tau^{i+1}, \nabla \varphi)_{K} \stackrel{(5.12\text{d})}{\leq} \eta_{\text{quad,F},K}^{i} \|(a^{i})^{\frac{1}{2}} \nabla \varphi\|_{K} \stackrel{(8.5)}{\leq} \eta_{\text{quad,F},K}^{i} \|\varphi\|_{i,K}. \quad (8.12) $$

Denote $t^{i} := (I - \Pi_{\ell,p})(S^{i} + L^{i} \phi_{\ell}^{i+1})$. To estimate the fourth term of (8.11), there are two possibilities. First, using the Cauchy–Schwarz inequality, (4.5a), and (8.5),

$$ (t^i, \varphi)_K \le \|(L_m^i|K)^{-\frac{1}{2}} t^i\|_K \|(L_m^i|K)^{\frac{1}{2}} \varphi\|_K \le (L_m^i|K)^{-\frac{1}{2}} \|t^i\|_K \|\varphi\|_{i,K}. \quad (8.13\text{a}) $$

Second, since the definition of $\Pi_{\ell,p}$ from (8.1a) allows to subtract a constant on each $K \in \mathcal{T}_{\ell}$, using again the Cauchy–Schwarz inequality, (4.5b), and (8.5)

$$ \begin{aligned} (t^i, \varphi)_K &= (t^i, \varphi - \frac{1}{|K|} \int_K \varphi)_K \stackrel{(2.1)}{\le} \frac{h_K}{\pi} \|t^i\|_K \|\nabla \varphi\|_K \\ &\le \frac{h_K}{\pi (a_m^i|K)^{\frac{1}{2}}} \|t^i\|_K \|(a_m^i|K)^{\frac{1}{2}} \nabla \varphi\|_K \le \frac{h_K}{\pi (a_m^i|K)^{\frac{1}{2}}} \|t^i\|_K \|\varphi\|_{i,K}. \end{aligned} \quad (8.13\text{b}) $$

Overall, using the notation $\tilde{w}_K^i$ from (8.3b),

$$ (t^i, \varphi)_K \le \tilde{w}_K^i \|t^i\|_K \|\varphi\|_{i,K} \stackrel{(5.12c)}{=} \eta_{\text{osc},K}^i \|\varphi\|_{i,K}. \quad (8.14) $$

The first term of (8.11) is simply estimated as (note that if $L^i$ is zero, this term vanishes)

$$ (L^i(u_{\ell}^{i+1} - \phi_{\ell}^{i+1}), \varphi)_K \le \|(L^i)^{\frac{1}{2}}(u_{\ell}^{i+1} - \phi_{\ell}^{i+1})\|_K \|(L^i)^{\frac{1}{2}}\varphi\|_K \stackrel{(5.12\text{b}), (8.5)}{\le} \eta_{\text{S},K}^i \|\varphi\|_{i,K}. \quad (8.15) $$

For the second term of (8.11), one can estimate as above

$$ (\Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + \sigma_{\ell}^{i+1}, \nabla \varphi)_{K} \le \|(a_{\text{m}}^{i})^{-\frac{1}{2}} (\Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + \sigma_{\ell}^{i+1})\|_{K} \|(a_{\text{m}}^{i})^{\frac{1}{2}} \nabla \varphi\|_{K}. \quad (8.16\text{a}) $$

However, this is not an optimal bound if reaction is the dominant effect (in a singularly perturbed regime which occurs, e.g., when $\tau \ll 1$ in Table 2). Nevertheless, since $(\Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + \sigma_{\ell}^{i+1})|_{K}$
is a polynomial lying in $RT_p(K)$ for all mesh elements $K$, it is possible to get using the trace and inverse inequalities (8.2) the bound

$$ (\Pi_{\ell,p-1}^{RT} \tau^{i+1} + \sigma_{\ell}^{i+1}, \nabla \varphi)_K \leq \eta_{F,K}^{i} \|\varphi\|_{i,K}. \quad (8.16b) $$

The proof of this estimate is given in [60, proof of Theorem 3.1], with $(a_m^i)^{1/2}$ replaced by $\epsilon$ and $(L_m^i)^{1/2}$ by $\kappa$, yielding the cut-off $w^i|_K$ given by (8.3a).

Combining (8.12)-(8.16) in (8.11) yields

$$ \langle \mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1}), \varphi \rangle \leq \sum_{K \in \mathcal{T}_{\ell}} (\eta_{F,K}^{i} + \eta_{S,K}^{i} + \eta_{\text{osc},K}^{i} + \eta_{\text{quad,F},K}^{i}) \|\varphi\|_{i,K}. \quad (8.17) $$

Upon using the Cauchy-Schwarz inequality one obtains (8.10) after noting that $\|\varphi\|_{1,u_{\ell}^i} = (\sum_{K \in \mathcal{T}_{\ell}} \|\varphi\|_{i,K}^2)^{1/2}$. $\square$

## 8.3 Local efficiency bounds

The efficiency proofs follow the methodology of [60]. We first prove the following lemma:

**Lemma 8.1 (Local efficiency bounds).** Under the assumptions of Theorem 5.7, one has

$$ [\eta_{\text{disc},\omega}^{i}]^2 \lesssim \sum_{\substack{\mathbf{a} \in \mathcal{V}_{\ell}, \\ \omega_{\mathbf{a}} \subseteq \tilde{\omega}}} \left[ \vartheta_{\mathbf{a},\ell,i}^2 \left( \sup_{\varphi \in H_0^1(\omega_{\mathbf{a}})} \frac{\langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle}{\|\varphi\|_{1,u_{\ell}^i,\omega_{\mathbf{a}}}} \right)^2 + [\eta_{\text{osc},\omega_{\mathbf{a}}}^i]^2 + \vartheta_{\mathbf{a},\ell,i}^2 [\eta_{\text{quad,F},\omega_{\mathbf{a}}}^i]^2 + [\eta_{\text{quad,S},\omega_{\mathbf{a}}}^i]^2 + [\eta_{\text{quad,S,osc},\omega_{\mathbf{a}}}^i]^2 \right], $$

where the constant in $\lesssim$ depends only on the space dimension $d$, the shape-regularity constant of $\mathcal{T}_\ell$, and the polynomial degree $p$.

*Proof.* Let $\mathcal{V}_{\tilde{\omega}} \subseteq \mathcal{V}_{\ell}$ be the collection of vertices of $\mathbf{a} \in \mathcal{V}_{\ell}$ such that $\omega_{\mathbf{a}} \subseteq \tilde{\omega}$. Then, one has

$$ \begin{aligned} [\eta_{\text{disc},\omega}^{i}]^2 &\stackrel{(5.13a)}{\lesssim} \sum_{K \in \mathcal{T}_{\ell}, K \subseteq \omega} ([\eta_{\text{osc},K}^{i}]^2 + [\eta_{\text{quad,F},K}^{i}]^2) + \sum_{K \in \mathcal{T}_{\ell}, K \subseteq \omega} ([\eta_{S,K}^{i}]^2 + [\eta_{F,K}^{i}]^2) \\ &\stackrel{(5.12)}{\lesssim} \sum_{\mathbf{a} \in \mathcal{V}_{\tilde{\omega}}} ([\eta_{\text{osc},\omega_{\mathbf{a}}}^{i}]^2 + [\eta_{\text{quad,F},\omega_{\mathbf{a}}}^{i}]^2) + \sum_{K \in \mathcal{T}_{\ell}, K \subseteq \omega} \left[ \|(L^i)^{\frac{1}{2}}(u_{\ell}^{i+1} - \phi_{\ell}^{i+1})\|_{K}^{2} \right. \\ &\qquad \left. + \|w^i(a_m^i)^{-\frac{1}{2}}(\Pi_{\ell,p-1}^{RT} \tau^{i+1} + \sigma_{\ell}^{i+1})\|_{K}^{2} \right]. \end{aligned} \quad (8.18a) $$

Then, calling the last summation above $T_1$, using the linearity of $\Pi_{\ell,p}$ from (8.1a), we have

$$ \begin{aligned} T_1 &\stackrel{(5.8),(8.9)}{=} \sum_{K \in \mathcal{T}_{\ell}, K \subseteq \omega} \left[ \left\| w^i(a_m^i)^{-\frac{1}{2}} \sum_{\mathbf{a} \in \mathcal{V}_K} (\psi_{\mathbf{a}} \Pi_{\ell,p-1}^{RT} \tau^{i+1} + \sigma_{\mathbf{a}}^{i+1}) \right\|_K^2 \right. \\ &\qquad \left. + \left\| (L^i)^{\frac{1}{2}} \sum_{\mathbf{a} \in \mathcal{V}_K} (\Pi_{\ell,p}(\psi_{\mathbf{a}} u_{\ell}^{i+1}) - \phi_{\mathbf{a}}^{i+1}) \right\|_K^2 \right] \\ &\lesssim \sum_{\mathbf{a} \in \mathcal{V}_{\tilde{\omega}}} \left[ \|w^i(a_m^i)^{-\frac{1}{2}}(\psi_{\mathbf{a}} \Pi_{\ell,p-1}^{RT} \tau^{i+1} + \sigma_{\mathbf{a}}^{i+1})\|_{\omega_{\mathbf{a}}}^2 \right. \\ &\qquad \left. + \|(L^i)^{\frac{1}{2}}(\Pi_{\ell,p}(\psi_{\mathbf{a}} u_{\ell}^{i+1}) - \phi_{\mathbf{a}}^{i+1})\|_{\omega_{\mathbf{a}}}^2 \right]. \end{aligned} \quad (8.18b) $$

The terms inside the summation now correspond to the minimization (5.7) and can be estimated as in the proof of [60, Theorem 4.4] by considering two cases. Let $\mathbf{a} \in \mathcal{V}_{\tilde{\omega}}$ be fixed. We henceforth
suppose that $L^i \neq 0$ on $\omega_a$, so that $L_M^i|_{\omega_a} > 0$ and $Q_a = P_p(T_a)$ from (5.6a); the pure diffusion case $L^i = 0$ on $\omega_a$ is standard (easier) and treated as in [24, 25] and the references therein.

**Diffusion-dominated case** ($h_{\omega_a}^2 \max_{\omega_a} L_M^i / \max_{\omega_a} a_M^i \le 1$). We adjust the methodology from [60], proof of Lemma 4.2, Case 2. Let $(v_a, q_a) \in V_a \times Q_a$ be given by

$$ q_a := (\Pi_{\ell,p,L^i}(\psi_a u_\ell^{i+1}))|_{\omega_a}, \\ v_a := \underset{\substack{\psi_\ell \in V_a, \\ \nabla \cdot w_\ell = - \Pi_{\ell,p}(\psi_a(S^i + L^i u_\ell^{i+1})) - \nabla \psi_a \cdot \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} \\ =: g_a}}{\operatorname{arg\,min}} \|\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + w_\ell\|_{\omega_a}. $$

The above minimization corresponds to a pure diffusion problem and has piecewise polynomial data $\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} \in \mathcal{RT}_p(T_a)$ and $g_a \in P_p(T_a)$. Moreover, for any interior vertex $a \in V_\ell^{\text{int}}$, we see from the second and third equalities in (8.8) that the mean value of the above divergence constraint is zero, $(g_a, 1)_{\omega_a} = 0$ so that the Neumann compatibility condition for the above $v_a \in V_a$ is satisfied (recall the zero normal trace constraint from (5.6b)). Using

$$ \Pi_{\ell,p}(L^i q_a) = \Pi_{\ell,p}(L^i \psi_a u_\ell^{i+1}), \quad (8.19) $$

the pair $(v_a, q_a)$ crucially satisfies the constraint in (5.7) in that

$$ \nabla \cdot v_a + \Pi_{\ell,p}(L^i q_a) = -\Pi_{\ell,p}(\psi_a S^i) - \nabla \psi_a \cdot \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}. \quad (8.20) $$

To see (8.19), note that from (8.1c), $q_a \in Q_a$ can be equivalently defined by

$$ (L^i q_a, q_\ell)_{\omega_a} = (L^i \psi_a u_\ell^{i+1}, q_\ell)_{\omega_a} \quad \forall q_\ell \in Q_a. $$

In particular

$$ q_a = \Pi_{\ell,p}(\psi_a u_\ell^{i+1}) $$

when $L^i$ is elementwise constant, cf. the discussion below (8.1c).

Define the quadrature contribution

$$ \eta_{\text{quad,S,a}}^i := \|(L^i)^{\frac{1}{2}}(\Pi_{\ell,p}(\psi_a u_\ell^{i+1}) - \Pi_{\ell,p,L^i}(\psi_a u_\ell^{i+1}))\|_{\omega_a} $$

and note that it vanishes when $L^i$ is elementwise constant. Let

$$ H_*^1(\omega_a) := \begin{cases} \{v \in H^1(\omega_a), (v,1)_{\omega_a} = 0\} & \text{if } a \in V_\ell^{\text{int}}, \\ \{v \in H^1(\omega_a), v = 0 \text{ on } \partial\omega_a \cap \{\psi_a > 0\}\} & \text{if } a \in V_\ell^{\text{ext}}. \end{cases} $$

Bounding from above the patchwise contribution of $T_1$ in (8.18b), since $(\sigma_a^{i+1}, \phi_a^{i+1})$ from (5.7) are the minimizers, we have, also relying on $w^i \le 1$ from (8.3a) ($w^i = 1$ when diffusion dominates),

$$ \begin{aligned} & \|w^i(a_m^i)^{-\frac{1}{2}}(\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + \sigma_a^{i+1})\|_{\omega_a}^2 + \|(L^i)^{\frac{1}{2}}(\Pi_{\ell,p}(\psi_a u_\ell^{i+1}) - \phi_a^{i+1})\|_{\omega_a}^2 \\ & \le \|w^i(a_m^i)^{-\frac{1}{2}}(\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + v_a)\|_{\omega_a}^2 + [\eta_{\text{quad,S,a}}^i]^2 \\ & \le \frac{1}{\min_{\omega_a} a_m^i} \|\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + v_a\|_{\omega_a}^2 + [\eta_{\text{quad,S,a}}^i]^2 \\ & \lesssim \frac{1}{\min_{\omega_a} a_m^i} \left[ \sup_{\varphi \in H_*^1(\omega_a)} \frac{(g_a, \varphi)_{\omega_a} - (\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, \nabla \varphi)_{\omega_a}}{\|\nabla \varphi\|_{\omega_a}} \right]^2 + [\eta_{\text{quad,S,a}}^i]^2 \\ & = \frac{1}{\min_{\omega_a} a_m^i} \left[ \sup_{\varphi \in H_*^1(\omega_a)} \frac{(\Pi_{\ell,p}(\psi_a(S^i + L^i u_\ell^{i+1})), \varphi)_{\omega_a} + (\Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1}, \nabla(\psi_a \varphi))_{\omega_a}}{\|\nabla \varphi\|_{\omega_a}} \right]^2 + [\eta_{\text{quad,S,a}}^i]^2 \\ & \le \frac{1}{\min_{\omega_a} a_m^i} \left[ \sup_{\varphi \in H_*^1(\omega_a)} \left( \frac{(S^i + L^i u_\ell^{i+1}, \psi_a \varphi)_{\omega_a} + (\tau^{i+1}, \nabla(\psi_a \varphi))_{\omega_a}}{\|\psi_a \varphi\|_{1, u_\ell^i, \omega_a}} \frac{\|\psi_a \varphi\|_{1, u_\ell^i, \omega_a}}{\|\nabla \varphi\|_{\omega_a}} \right) \right. \\ & \qquad \left. + \sup_{\varphi \in H_*^1(\omega_a)} \frac{((\Pi_{\ell,p} - I)(\psi_a(S^i + L^i u_\ell^{i+1})), \varphi)_{\omega_a} + ((\Pi_{\ell,p-1}^{\text{RT}} - I)\tau^{i+1}, \nabla(\psi_a \varphi))_{\omega_a}}{\|\nabla \varphi\|_{\omega_a}} \right]^2 + [\eta_{\text{quad,S,a}}^i]^2. \end{aligned} \quad (8.21a) $$
The third (crucial) inequality follows from [60, Lemma 4.1], itself evoking [9, Theorem 7] and [25, Corollaries 3.3 and 3.6]. In the last expression in (8.21a), recalling (5.2) and (8.4), the first term, say $T_2$, features the residual $\langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \psi_a \varphi \rangle$; note also that $\psi_a \varphi \in H_0^1(\Omega)$ from the definition of the space $H_*^1(\omega_a)$. The second term $T_3$ is then related to quadrature and oscillation errors. We now bound them separately.

We first bound $T_2$. For this purpose, note that, cf., e.g., [9, estimate below (11)],

$$ \|\nabla(\psi_a\varphi)\|_{\omega_a} \lesssim \|\nabla\varphi\|_{\omega_a} \quad \forall \varphi \in H_*^1(\omega_a). $$

Moreover, the Poincaré inequality (2.1) and (4.5) give

$$ \begin{aligned} \|\psi_a \varphi\|_{1, u_{\ell}^i, \omega_a}^2 &= \|(L^i)^{\frac{1}{2}}(\psi_a \varphi)\|_{\omega_a}^2 + \|(a^i)^{\frac{1}{2}} \nabla(\psi_a \varphi)\|_{\omega_a}^2 \\ &\lesssim (h_{\omega_a}^2 \max_{\omega_a} L_M^i + \max_{\omega_a} a_M^i) \|\nabla(\psi_a \varphi)\|_{\omega_a}^2. \end{aligned} $$

Now, straightforwardly, under the assumption $h_{\omega_a}^2 \max_{\omega_a} L_M^i / \max_{\omega_a} a_M^i \le 1$ of (5.15a),

$$ \frac{1}{\min_{\omega_a} a_m^i} (h_{\omega_a}^2 \max_{\omega_a} L_M^i + \max_{\omega_a} a_M^i) \le \left( h_{\omega_a}^2 \frac{\max_{\omega_a} L_M^i}{\max_{\omega_a} a_M^i} + 1 \right) \frac{\max_{\omega_a} a_M^i}{\min_{\omega_a} a_m^i} \le 2 \frac{\max_{\omega_a} a_M^i}{\min_{\omega_a} a_m^i}. $$

Consequently,

$$ \begin{aligned} T_2 &= \frac{1}{\min_{\omega_a}(a_m^i)^{\frac{1}{2}}} \sup_{\varphi \in H_0^1(\omega_a)} \left( \frac{\langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \psi_a \varphi \rangle}{\|\psi_a \varphi\|_{1, u_{\ell}^i, \omega_a}} \frac{\|\psi_a \varphi\|_{1, u_{\ell}^i, \omega_a}}{\|\nabla \varphi\|_{\omega_a}} \right) \\ &\lesssim \vartheta_{a, \ell, i} \sup_{\varphi \in H_0^1(\omega_a)} \frac{\langle \mathcal{R}_{\text{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^i, \omega_a}}. \end{aligned} \quad (8.21b) $$

As for $T_3$, note that, using the above bound $\|\nabla(\psi_a \varphi)\|_{\omega_a} \lesssim \|\nabla\varphi\|_{\omega_a}$,

$$ \begin{aligned} & \frac{1}{\min_{\omega_a}(a_m^i)^{\frac{1}{2}}} \left( \left( \Pi_{\ell,p-1}^{\text{RT}} - \mathbb{I} \right) \tau^{i+1}, \nabla(\psi_a \varphi) \right)_{\omega_a} \\ & \stackrel{(5.15a)}{\le} \vartheta_{a, \ell, i} \|(a^i)^{-\frac{1}{2}} (\Pi_{\ell,p-1}^{\text{RT}} - \mathbb{I}) \tau^{i+1} \|_{\omega_a} \|\nabla(\psi_a \varphi)\|_{\omega_a} \stackrel{(5.12d)}{\lesssim} \vartheta_{a, \ell, i} \eta_{\text{quad,F},\omega_a}^i \|\nabla\varphi\|_{\omega_a}. \end{aligned} $$

Let

$$ \eta_{\text{quad,S,osc},a}^i := \frac{h_{\omega_a}}{\min_{\omega_a}(a_m^i)^{\frac{1}{2}}} \|(\Pi_{\ell,p} - I)(\psi_a(S^i + L^i u_{\ell}^{i+1}))\|_{\omega_a}. $$

Since $((\Pi_{\ell,p} - I)\zeta, \varphi)_{\omega_a} = ((\Pi_{\ell,p} - I)\zeta, \varphi - \Pi_{\ell,p}\varphi)_{\omega_a}$ for $\zeta \in L^2(\omega_a)$, the Poincaré inequality (2.1) gives

$$ \frac{1}{\min_{\omega_a}(a_m^i)^{\frac{1}{2}}} \left( \left( \Pi_{\ell,p} - I \right) \left( \psi_a(S^i + L^i u_{\ell}^{i+1}) \right), \varphi \right)_{\omega_a} \stackrel{(2.1)}{\lesssim} \eta_{\text{quad,S,osc},a}^i \|\nabla\varphi\|_{\omega_a}. $$

Thus, overall,

$$ T_3 \lesssim \vartheta_{a, \ell, i} \eta_{\text{quad,F}, \omega_a}^i + \eta_{\text{quad,S,osc},a}^i. \quad (8.21c) $$

The claim of Lemma 8.1 in this case follows by combining the inequalities in (8.18) and (8.21).

**Reaction-dominated case** ($h_{\omega_a}^2 \max_{\omega_a} L_M^i / \max_{\omega_a} a_M^i > 1$). We adjust the methodology from [60], proof of Lemma 4.2, Case 1. Hence, one chooses $(v_a, q_a) \in V_a \times Q_a$ as

$$ (L^i q_a, q_\ell)_{\omega_a} = (-\Pi_{\ell,p}(\psi_a S^i) + \psi_a \nabla \cdot \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} - \rho_a, q_\ell)_{\omega_a} \quad \forall q_\ell \in Q_a, \quad (8.22a) $$

$$ v_a := \arg \min_{\substack{w_\ell \in V_a, \\ \nabla \cdot w_\ell = g_a}} \|(\psi_a \Pi_{\ell,p-1}^{\text{RT}} \tau^{i+1} + w_\ell)\|_{\omega_a}, \quad (8.22b) $$
where $\nabla \cdot_{\ell}$ is the elementwise divergence and

$$ \rho_a := \frac{1}{|\omega_a|} (-S^i - L^i u_{\ell}^{i+1} + \nabla \cdot_{\ell} \Pi_{\ell, p-1}^{RT} \tau^{i+1}, \psi_a)_{\omega_a}, $$

$$ g_a := -\nabla \cdot_{\ell} (\psi_a \Pi_{\ell, p-1}^{RT} \tau^{i+1}) + \rho_a. $$

Observe that, $(g_a, 1)_{\omega_a} = 0$ using (4.9); we also need (8.1b) and $\nabla \psi_a \in \mathcal{RT}_0(\mathcal{T}_{\ell}; \mathbb{R}^d)$. As in (8.20), we also have

$$ \nabla \cdot v_a + \Pi_{\ell,p}(L^i q_a) = -\Pi_{\ell}(\psi_a S^i) - \nabla \psi_a \cdot \Pi_{\ell, p-1}^{RT} \tau^{i+1}. $$

Hence, an upper bound of the patch-wise terms in (8.18b) can be obtained from Definition 5.4 with this particular choice of $(v_a, q_a)$. The proof follows [60, proof of Lemma 4.2, Case 1], itself relying on [64]. □

## 8.4 Proof of efficiency estimates in Theorem 5.7

We can now accomplish the proof of the efficiency bounds in Theorem 5.7. Observe that $\vartheta_{a, \ell, i}^2 \le \vartheta_{\tilde{\omega}, \ell, i}^2$ by definition and

$$ \sum_{\substack{a \in \mathcal{V}_{\ell, i} \\ \omega_a \subseteq \tilde{\omega}}} ( [\eta_{\mathrm{osc}, \omega_a}^i]^2 + \vartheta_{a, \ell, i}^2 [\eta_{\mathrm{quad}, F, \omega_a}^i]^2 + [\eta_{\mathrm{quad}, S, \omega_a}^i]^2 + [\eta_{\mathrm{quad}, S, \mathrm{osc}, \omega_a}^i]^2 ) \\ \lesssim [\eta_{\mathrm{osc}, \tilde{\omega}}^i]^2 + \vartheta_{\tilde{\omega}, \ell, i}^2 [\eta_{\mathrm{quad}, F, \tilde{\omega}}^i]^2 + [\eta_{\mathrm{quad}, S, \tilde{\omega}}^i]^2 + [\eta_{\mathrm{quad}, S, \mathrm{osc}, \tilde{\omega}}^i]^2. $$

Hence, the proof of (5.14b) is completed from Lemma 8.1 by showing that

$$ \sum_{\substack{a \in \mathcal{V}_{\ell} \\ \omega_a \subseteq \tilde{\omega}}} \left( \sup_{\varphi \in H_0^1(\omega_a)} \frac{\langle \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^i, \omega_a}} \right)^2 \lesssim \| \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}) \|_{-1, u_{\ell}^i, \tilde{\omega}}^2 \quad (8.23) $$

followed by adding $[\eta_{\mathrm{lin}, \omega}^i]^2$ to both sides of the inequality (see the definition of $\eta_{\tilde{\omega}}^i$ in (5.13c)).

To show (8.23), we proceed as in [16, Theorem 3.3] and the references therein. For each $a \in \mathcal{V}_{\ell}$ such that $\omega_a \subseteq \tilde{\omega}$, let us consider the Riesz representer $r_a \in H_0^1(\omega_a) \subseteq H_0^1(\tilde{\omega})$ satisfying for all $\varphi \in H_0^1(\omega_a)$,

$$ ((r_a, \varphi))_{u_{\ell}^i} = \langle \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle, \quad \text{so that} \quad \sup_{\varphi \in H_0^1(\omega_a)} \frac{\langle \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^i, \omega_a}} = \|r_a\|_{1, u_{\ell}^i, \omega_a}. \quad (8.24) $$

Then, defining $r := \sum_{a \in \mathcal{V}_{\ell}, \omega_a \subseteq \tilde{\omega}} r_a \in H_0^1(\tilde{\omega})$, one has

$$ \sum_{\substack{a \in \mathcal{V}_{\ell} \\ \omega_a \subseteq \tilde{\omega}}} \|r_a\|_{1, u_{\ell}^i, \omega_a}^2 = \sum_{\substack{a \in \mathcal{V}_{\ell} \\ \omega_a \subseteq \tilde{\omega}}} \langle \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), r_a \rangle = \langle \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}), r \rangle \\ \le \| \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}) \|_{-1, u_{\ell}^i, \tilde{\omega}} \|r\|_{1, u_{\ell}^i, \tilde{\omega}} \le (d+1)^{\frac{1}{2}} \| \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}) \|_{-1, u_{\ell}^i, \tilde{\omega}} \left( \sum_{\substack{a \in \mathcal{V}_{\ell} \\ \omega_a \subseteq \tilde{\omega}}} \|r_a\|_{1, u_{\ell}^i, \omega_a}^2 \right)^{\frac{1}{2}}, $$

Finally, the global efficiency estimate (5.14c) follows from choosing $\omega = \Omega$ in (5.14b) together with $\| \mathcal{R}_{\mathrm{disc}}^{u_{\ell}^i}(u_{\ell}^{i+1}) \|_{-1, u_{\ell}^i}^2 + [\eta_{\mathrm{lin}, \Omega}^i]^2 = \| \mathcal{R}(u_{\ell}^{i+1}) \|_{-1, u_{\ell}^i}^2$ from Theorem 5.3 and $\vartheta_{\Omega, \ell, i} \ge 1$ from (5.15).
# Appendices

## Appendix A Iterative linearizations in abstract spaces: well-posedness and consistency

Here, we generalize the iterative linearization defined in Sections 4.3 and 4.4 to an abstract setting, and show well-posedness and consistency of the iterations. With reference to the classes of problems discussed in Section 3, let us introduce the space

$$W_q := \begin{cases} W_0^{1,q}(\Omega) & \text{for gradient-dependent diffusivity case (3.1),} \\ L^\infty(\Omega) \cap W_0^{1,q}(\Omega) & \text{for gradient-independent diffusivity case (3.3),} \end{cases} \qquad (A.1)$$

for $q \in [2, \infty]$. Observe that the finite element space $V_\ell \subset W_q$ due to (4.2). We equip $W_q$ with the norm $\|v\|_{W_q} := \|\nabla v\|_{L^q(\Omega)}$ for (3.1), and $\|v\|_{W_q} := \|\nabla v\|_{L^q(\Omega)} + \|v\|_{L^\infty(\Omega)}$ for (3.3).

In extension of (4.4), for a given parameter function $v \in W_2$, let

$$((\xi, \zeta))_v := (L(\cdot, v)\xi, \zeta) + (a(\cdot, v)\nabla\xi, \nabla\zeta), \quad \xi, \zeta \in H_0^1(\Omega). \qquad (A.2)$$

Similarly, extending (4.15) and (4.18), let

$$\|\xi\|_{1,v} := ((\xi, \xi))_v^{\frac{1}{2}}, \quad \xi \in H_0^1(\Omega), \qquad (A.3a)$$

$$\|\zeta\|_{-1,v} := \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \zeta, \varphi \rangle}{\|\varphi\|_{1,v}}, \quad \zeta \in H^{-1}(\Omega). \qquad (A.3b)$$

Generalizing Assumption 4.1, here we assume

**Assumption A.1 (Coefficients L and a).** The function $L: \Omega \times \mathbb{R} \to \mathbb{R}$ and the symmetric tensor $a: \Omega \times \mathbb{R} \to \mathbb{R}^d \times \mathbb{R}^d$ are measurable and satisfy for continuous functions $L_m, L_M: [0, \infty) \to [0, \infty)$ and $a_m, a_M: [0, \infty) \to (0, \infty)$, for almost every $x \in \Omega$, and for all $v \in W_2$,

$$0 \le L_m(\|v\|_{W_2}) \le L(x, v(x)) \le L_M(\|v\|_{W_2}), \qquad (A.4a)$$

$$0 < a_m(\|v\|_{W_2})|y|^2 \le y^T a(x, v(x)) y \le a_M(\|v\|_{W_2})|y|^2 \quad \forall y \in \mathbb{R}^d \setminus \{0\}. \qquad (A.4b)$$

The functions $L_{m/M}$ and $a_{m/M}$ can be completely independent of the argument, for instance for the Zarantonello scheme of Table 1. Here, we adhere to a more broad setting of (A.4).

We show that the linearization scheme (4.7) is a special case of a more general iterative linearization strategy involving the inner product $((\cdot, \cdot))_v$ from (A.2), and we also show that it is a well-posed and consistent strategy. For this purpose, let us introduce

**Definition A.2 (Linearization operator $\mathcal{I}$).** Let $((\cdot, \cdot))_v$ be defined in (A.2), with the coefficients $L$ and $a$ satisfying Assumption A.1. Let $q \ge 2$ be fixed. Then, we define the linearization operator $\mathcal{I}: W_q \to W_q$ as: for a given $v \in W_q$, $\mathcal{I}v \in W_q$ solves

$$((\mathcal{I}v - v, \varphi))_v = -\langle \mathcal{R}(v), \varphi \rangle \quad \forall \varphi \in H_0^1(\Omega). \qquad (A.5)$$

Here $\mathcal{I}v - v$ is the increment in $v$ towards the weak solution $u \in H_0^1(\Omega)$ of (1.1), and $\langle \mathcal{R}(v), \varphi \rangle$ defines the residual relation with respect to $v$. With reference to (4.7) and (4.13), $u_{\langle \ell \rangle}^{i+1} = \mathcal{I}u_\ell^i$. The operator $\mathcal{I}$ extends the linear iteration in a continuous level for some initial guess $u_0 \in W_q$, i.e. $u^i = (\mathcal{I})^i u_0$ for $i \ge 1$. We need to additionally assume here that the space dimension is limited to $1 \le d \le 3$ to ensure its well-posedness:
**Proposition A.3** (Well-posedness of the operator $\mathcal{I}$). Assume (A1)-(A2) if $\mathcal{R}: H_0^1(\Omega) \to H^{-1}(\Omega)$ is defined by (3.1), and (B1)-(B5) if $\mathcal{R}$ is defined by (3.3). For spatial dimension $d \in \{1, 2, 3\}$, let a finite $q \ge 2$ satisfy $d < q \le 2 \max(2, d)/(\max(2, d) - 2)$. Let the space $W_q$ be defined in (A.1), and let $((\cdot, \cdot))_v$ be defined in (A.2), with the coefficients $L$ and $\mathbf{a}$ satisfying Assumption A.1. Then, the linearization operator $\mathcal{I}: W_q \to W_q$ of Definition A.2 is well-defined.

*Proof.* Let $v \in W_q$ be fixed. Observe that this implies $v \in W_2$ since $W_q \subseteq W_2$ as $q \ge 2$, and in particular $v \in H_0^1(\Omega)$. Consequently, in view of Assumption A.1, $((\cdot, \cdot))_v$ is a $H_0^1(\Omega)$-inner product. Then, existence and uniqueness of $\mathcal{I}v \in H_0^1(\Omega)$ follow from the Riesz representation theorem, and we are left to show that $\mathcal{I}v \in W_q$.

First, let $q = 2$. By the assumption $d < q$, this only arises when $d = 1$. In this case, cf., e.g., [26, Chapter 5.6], $H_0^1(\Omega) \subset C(\bar{\Omega}) \subset L^\infty(\Omega)$, so that indeed $\mathcal{I}v \in W_q$.

Second, let $q > \max\{2, d\}$ be finite and consider the expression (3.1) of $\mathcal{R}$. Since $\mathcal{I}v \in H_0^1(\Omega)$ and $q \le 2 \max(2, d)/(\max(2, d) - 2)$, Sobolev's embedding theorem [26, Chapter 5.6] further gives $\mathcal{I}v \in L^q(\Omega)$. Using (A.4a), we then have $L(\mathbf{x}, v) \mathcal{I}v \in L^q(\Omega)$. Moreover, one has from Assumption 3.1 that

$$ \begin{aligned} |\sigma(\mathbf{x}, \nabla v)| &\le |\sigma(\mathbf{x}, \nabla v) - \sigma(\mathbf{x}, 0)| + |\sigma(\mathbf{x}, 0)| &\stackrel{(A1)}{\le} \sigma_M |\nabla v| + |\sigma(\mathbf{x}, 0)| &\stackrel{(A1)}{\in} L^q(\Omega), \\ |f(\mathbf{x}, v)| &\le |f(\mathbf{x}, v) - f(\mathbf{x}, 0)| + |f(\mathbf{x}, 0)| &\stackrel{(A2)}{\le} f_M |v| + |f(\mathbf{x}, 0)| &\stackrel{(A2)}{\in} L^q(\Omega). \end{aligned} $$

Thus, rearranging (A.5) gives

$$ (\mathbf{a}(\cdot, v)\nabla \mathcal{I}v, \nabla \varphi) = (\mathcal{S}^*, \varphi) + (\mathbf{F}^*, \nabla \varphi) \quad \forall \varphi \in H_0^1(\Omega), \text{ with } \mathcal{S}^* \in L^q(\Omega) \text{ and } \mathbf{F}^* \in L^q(\Omega; \mathbb{R}^d). $$

Now, as $\mathbf{a}(\mathbf{x}, v)$ satisfies the ellipticity condition due to (A.4b), applying Meyers' theorem [51, Theorem 1], we see that $\mathcal{I}v \in W_0^{1,q}(\Omega) = W_q$.

Third, let $q > \max\{2, d\}$ be finite and consider the expression (3.3) of $\mathcal{R}$. Using Assumption 3.3, this yields

$$ f(\cdot, v) \in L^q(\Omega) \text{ and } \bar{\mathbf{K}}(\cdot)(\mathcal{D}(\cdot, v)\nabla v + \mathbf{q}(\cdot, v)) \in L^q(\Omega; \mathbb{R}^d). $$

The last inclusion holds since $\bar{\mathbf{K}}\mathbf{q} \in L^\infty(\Omega; \mathbb{R}^d)$, and $\mathcal{D}\bar{\mathbf{K}}\nabla v \in L^q(\Omega; \mathbb{R}^d)$ as $v \in W_0^{1,q}(\Omega)$. Then the conclusion $\mathcal{I}v \in W_0^{1,q}(\Omega)$ follows as above. Finally $\mathcal{I}v \in W^{1,q}(\Omega)$ and $q > d$ implies $\mathcal{I}v \in C(\bar{\Omega}) \subset L^\infty(\Omega)$ from the Morrey's inequality [26, Chapter 5.6], so that indeed $\mathcal{I}v \in W_q$. $\square$

**Proposition A.4** (Consistency of the operator $\mathcal{I}$). Under the assumptions of Proposition A.3, let $u \in W_2$ be the unique weak solution of (1.1) for a continuous operator $\mathcal{R}: H_0^1(\Omega) \to H^{-1}(\Omega)$. Then, $\mathcal{I}v = v$ if and only if $v = u$. Moreover, for a given $u_0 \in W_q$, if the sequence $\{(\mathcal{I})^i u_0 \in W_q\}$ is bounded uniformly in $W_2$ with respect to $i \ge 0$ and Cauchy in $H_0^1(\Omega)$, then $\text{dist}((\mathcal{I})^i u_0, u) \to 0$ as $i \to \infty$ with the metric dist: $H_0^1(\Omega) \times H_0^1(\Omega) \to [0, \infty)$ introduced in (1.2).

*Proof.* If $u \in W_2$ is the weak solution of (1.1), then $((\mathcal{I}u - u, \varphi))_u = -\langle \mathcal{R}(u), \varphi \rangle = 0$ for all $\varphi \in H_0^1(\Omega)$, yielding $\mathcal{I}u - u = 0$. On the other hand, if $\mathcal{I}v = v$ then also $\langle \mathcal{R}(v), \varphi \rangle = 0$ for all $\varphi \in H_0^1(\Omega)$ which gives $v = u$ due to the uniqueness of $u$.

If $(\mathcal{I})^i u_0$ is bounded in $W_2$, then there exists a constant $C > 0$ such that for all $\varphi \in H_0^1(\Omega)$,

$$ \begin{aligned} |\langle \mathcal{R}((\mathcal{I})^i u_0) - \mathcal{R}(u), \varphi \rangle| &= |\langle \mathcal{R}((\mathcal{I})^i u_0), \varphi \rangle| = |((\mathcal{I})^{i+1} u_0 - (\mathcal{I})^i u_0, \varphi))_{(\mathcal{I})^i u_0}| \\ &\le C \|\nabla ((\mathcal{I})^{i+1} u_0 - (\mathcal{I})^i u_0)\| \|\nabla \varphi\|. \end{aligned} $$

The right-hand side goes to 0 as $i \to \infty$ as the sequence $\{(\mathcal{I})^i u_0\}_{i \ge 0}$ is Cauchy. Consequently $\text{dist}((\mathcal{I})^i u_0, u) \to 0$ due to (1.2). $\square$

<figure><img src="image_15.png" /></figure>
## Appendix B Error estimates in a fixed norm

The norm $|||\cdot|||_{-1, u_{\ell}^i}$, introduced in (4.18), is not a fixed norm since its definition changes with the iteration index $i$. In this section, we show how the bounds for $|||\mathcal{R}(u_{\ell}^i)|||_{-1, u_{\ell}^i}$, presented in Section 5, can be converted to give estimates of $\mathcal{R}(u_{\ell}^i)$ in the fixed norm $|||\cdot|||_{-1, u_{\ell}}$, where $u_{\ell} \in V_{\ell}$ is the finite element solution of the nonlinear problem introduced in (4.3). For this purpose, we adopt the setting of Appendix A, in particular (A.2), (recall $V_{\ell} \subset W_q \subseteq W_2$ from (4.2) and since $q \ge 2$), and assume that:

**Assumption B.1** (Dependence of $((\cdot, \cdot))_v$ on $v \in V_{\ell}$). For $v, w \in V_{\ell}$ and $\xi \in H_0^1(\Omega)$, there exists a positive function $F \in C(\Omega \times \mathbb{R})$, Lipschitz with respect to the second argument $\xi \in \mathbb{R}$, i.e.,

$$|F(\mathbf{x}, \xi_1) - F(\mathbf{x}, \xi_2)| \le F_M |\xi_1 - \xi_2| \quad \forall \mathbf{x} \in \Omega \text{ and } \xi_1, \xi_2 \in \mathbb{R}, \quad (B.1)$$

such that

$$\|\xi\|_{1,w}^2 \le \max \left( \sup_{\mathbf{x} \in \Omega} \left\{ \frac{F(\mathbf{x},w)}{F(\mathbf{x},v)} \right\}, 1 \right) \|\xi\|_{1,v}^2. \quad (B.2)$$

**Remark B.2** (Applicability of Assumption B.1). Assumption B.1 is satisfied for most of the schemes presented in Tables 1 and 2 with the corresponding F-functions given in Table 3. For example, in the case of the L-scheme, one may choose $F = \mathcal{D}$ since

$$\begin{aligned} \|\xi\|_{1,w}^2 &= \|L^{\frac{1}{2}}\xi\|^2 + \|\tau^{\frac{1}{2}} \bar{K}^{\frac{1}{2}}(\cdot) D(\cdot, w)^{\frac{1}{2}} \nabla \xi\|^2 \\ &\le \max \left( \sup_{\mathbf{x} \in \Omega} \left\{ \frac{D(\mathbf{x},w)}{D(\mathbf{x},v)} \right\}, 1 \right) \left[ \|L^{\frac{1}{2}}\xi\|^2 + \|\tau^{\frac{1}{2}} \bar{K}^{\frac{1}{2}}(\cdot) D(\cdot, v)^{\frac{1}{2}} \nabla \xi\|^2 \right] \\ &= \max \left( \sup_{\mathbf{x} \in \Omega} \left\{ \frac{D(\mathbf{x},w)}{D(\mathbf{x},v)} \right\}, 1 \right) \|\xi\|_{1,v}^2. \end{aligned}$$

However, the Newton scheme and the Kačanov scheme do not satisfy Assumption B.1 directly. In these cases, the function $F$ also needs to depend on derivatives, i.e., $F = F(\mathbf{x}, v, \nabla v)$. We have not considered this option here for simplicity.

<table><caption>Table 3: The choices of $F$ in Assumption B.1 for different linearization schemes given in Tables 1 and 2. For Picard, Jäger-Kačur, and M-scheme, the choice yielding higher $\sup_{\mathbf{x} \in \Omega}\{F(\mathbf{x}, w)/F(\mathbf{x}, v)\}$ value needs to be selected.</caption><thead><tr><th>Scheme</th><th>$F(\mathbf{x}, \xi)$</th></tr></thead><tbody><tr><td>Zarantonello</td><td>$\Lambda$</td></tr><tr><td>Picard</td><td>$\partial_{\xi} f(\mathbf{x}, \xi)$ or $D(\mathbf{x}, \xi)$</td></tr><tr><td>Jäger-Kačur</td><td>$\max_{\zeta \in \mathbb{R}} \left( \frac{f(\mathbf{x}, \zeta) - f(\mathbf{x}, \xi)}{\zeta - \xi} \right)$ or $D(\mathbf{x}, \xi)$</td></tr><tr><td>L-scheme</td><td>$D(\mathbf{x}, \xi)$</td></tr><tr><td>M-scheme</td><td>$\partial_{\xi} f(\mathbf{x}, \xi) + M\tau$ or $D(\mathbf{x}, \xi)$</td></tr></tbody></table>

If stable iterative schemes converge with a linear rate in $H_0^1(\Omega)$, and more importantly in $L^2(\Omega)$, they are also expected to converge linearly in $L^\infty(\Omega)$. This is proven for the L- and the M-schemes in Proposition 3.2 of [52] for the continuous solutions. Here we assume it for the finite element solutions:

**Assumption B.3** (Contraction in $L^\infty$ of the linear iterations). There exists $\alpha \in (0, 1)$ such that

$$\|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{L^{\infty}(\Omega)} \le \alpha \|u_{\ell}^{i} - u_{\ell}^{i-1}\|_{L^{\infty}(\Omega)} \quad \forall i \ge 1. \quad (B.3)$$
**Remark B.4** (Estimating $\alpha \in (0,1)$). Note that the theoretical estimates for the contraction rates always depend on the Lipschitz/monotonicity ratio $\lambda_M/\lambda_m$ from (1.2), and thus give a pessimistic theoretical value of $\alpha \approx 1$ whenever this ratio is large. However, in Assumption B.3, the parameter $\alpha$ can be directly estimated from the numerical solutions as

$$ \alpha \ge \alpha_i := \|u_{\ell}^{i+1} - u_{\ell}^i\|_{L^{\infty}(\Omega)} / \|u_{\ell}^i - u_{\ell}^{i-1}\|_{L^{\infty}(\Omega)}. $$

Thus, for the guaranteed linearly converging schemes, like the Zarantonello, Jäger-Kačur, L- and M-schemes, one expects to get a tight practical bound on $\alpha \in (0,1)$ away from 1 with only a couple of iterations.

**Theorem B.5** (Error estimates in the fixed norm $\|\cdot\|_{-1, u_\ell}$). Let $u \in H_0^1(\Omega)$ solve (1.1) and $u_\ell \in V_\ell$ solve (4.3). Let the sequence $\{u_\ell^i\}_{i \ge 0} \subset V_\ell$ be defined by (4.7), using the scalar product (4.4) and the norms (4.15) and (4.18). Let Assumptions B.1 and B.3 hold. Then, for any $\tilde{u} \in H_0^1(\Omega)$,

$$ \frac{\|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}^i}}{(1 + \mathcal{C}_m^i \|u_{\ell}^{i+1} - u_{\ell}^i\|_{L^{\infty}(\Omega)})^{\frac{1}{2}}} \le \|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}} \le (1 + \mathcal{C}_m^i \|u_{\ell}^{i+1} - u_{\ell}^i\|_{L^{\infty}(\Omega)})^{\frac{1}{2}} \|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}^i}, \quad (B.4) $$

where the constants $\mathcal{C}_m^i, \mathcal{C}_M^i \ge 0$, given by

$$ \mathcal{C}_m^i := \frac{1}{1-\alpha} \sup_{\mathbf{x} \in \Omega} \left\{ \frac{1}{F(\mathbf{x}, u_{\ell}^i)} \sup_{v \in I_{\ell}^i(\mathbf{x})} \left| \frac{F(\mathbf{x}, v) - F(\mathbf{x}, u_{\ell}^i)}{v - u_{\ell}^i} \right| \right\}, \quad (B.5a) $$

$$ \mathcal{C}_M^i := \frac{1}{1-\alpha} \sup_{\mathbf{x} \in \Omega} \left\{ F(\mathbf{x}, u_{\ell}^i) \sup_{v \in I_{\ell}^i(\mathbf{x})} \left| \frac{\frac{1}{F(\mathbf{x},v)} - \frac{1}{F(\mathbf{x},u_{\ell}^i)}}{v - u_{\ell}^i} \right| \right\}, \quad (B.5b) $$

are uniformly bounded with respect to iteration $i$, where the interval $I_{\ell}^i(\mathbf{x})$ is defined for all $\mathbf{x} \in \Omega$ as

$$ I_{\ell}^{i}(\mathbf{x}) := \left[ u_{\ell}^{i}(\mathbf{x}) - \frac{\|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{L^{\infty}(\Omega)}}{1-\alpha}, u_{\ell}^{i}(\mathbf{x}) + \frac{\|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{L^{\infty}(\Omega)}}{1-\alpha} \right]. $$

*Proof.* Observe that Assumption B.3 implies that $u_{\ell}^i \to u_\ell$ as $i \to \infty$ and

$$ \|u_{\ell}^{n+1} - u_{\ell}^{n}\|_{L^{\infty}(\Omega)} \le \alpha^{n-i} \|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{L^{\infty}(\Omega)} $$

for $n \ge i$. This gives

$$ \begin{aligned} \|u_\ell - u_\ell^i\|_{L^\infty(\Omega)} &\le \sum_{n=i}^\infty \|u_\ell^{n+1} - u_\ell^n\|_{L^\infty(\Omega)} \le \|u_\ell^{i+1} - u_\ell^i\|_{L^\infty(\Omega)} \sum_{n=0}^\infty \alpha^n \\ &= \frac{\|u_\ell^{i+1} - u_\ell^i\|_{L^\infty(\Omega)}}{1-\alpha} \le \frac{\alpha^i}{1-\alpha} \|u_\ell^1 - u_\ell^0\|_{L^\infty(\Omega)}, \end{aligned} \quad (B.6a) $$

which also yields for all $\mathbf{x} \in \Omega$ that

$$ u_{\ell}^{i}(\mathbf{x}) - \frac{\|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{L^{\infty}(\Omega)}}{1-\alpha} \le u_{\ell}(\mathbf{x}) \le u_{\ell}^{i}(\mathbf{x}) + \frac{\|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{L^{\infty}(\Omega)}}{1-\alpha}, \quad \text{or} \quad u_{\ell}(\mathbf{x}) \in I_{\ell}^{i}(\mathbf{x}). \quad (B.6b) $$

Hence, if $\tilde{\varphi} = \arg\max_{\varphi \in H_0^1(\Omega)} (|\langle \mathcal{R}(\tilde{u}), \varphi \rangle| / \|\varphi\|_{1, u_\ell})$, then Assumption B.1 yields

$$ \|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}} = \frac{|\langle \mathcal{R}(\tilde{u}), \tilde{\varphi} \rangle|}{\|\tilde{\varphi}\|_{1, u_{\ell}}} = \left( \frac{|\langle \mathcal{R}(\tilde{u}), \tilde{\varphi} \rangle|}{\|\tilde{\varphi}\|_{1, u_{\ell}^i}} \right) \left( \frac{\|\tilde{\varphi}\|_{1, u_{\ell}^i}}{\|\tilde{\varphi}\|_{1, u_{\ell}}} \right) \le \|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}^i} \sqrt{\sup_{\Omega} \frac{F(u_{\ell}^i)}{F(u_{\ell})}}. \quad (B.7) $$

Here, we have dropped the $\mathbf{x}$-dependence of $F$ for simplifying the notation and we have assumed that $\sup_{\Omega} F(u_\ell^i)/F(u_\ell) > 1$, since the $\sup_{\Omega} F(u_\ell^i)/F(u_\ell) \le 1$ case is trivial in that it gives
$\|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}} \le \|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}^i}$ directly. Then, we get the upper estimate of $\|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}}$ by observing

$$
\begin{align}
\sup_{\Omega} \frac{F(u_{\ell}^i)}{F(u_{\ell})} &= 1 + \sup_{\Omega} \frac{F(u_{\ell}^i) - F(u_{\ell})}{F(u_{\ell})} = 1 + \sup_{\Omega} \left\{ F(u_{\ell}^i) \left( \frac{1}{F(u_{\ell})} - \frac{1}{F(u_{\ell}^i)} \right) \right\} \nonumber \\
&\stackrel{\text{(B.6b)}}{\le} 1 + \sup_{\Omega} \left\{ F(u_{\ell}^i) \left| \frac{1}{u_{\ell} - u_{\ell}^i} - \frac{1}{F(u_{\ell}^i)} \right| \right\} \|u_{\ell} - u_{\ell}^i\|_{L^{\infty}(\Omega)} \nonumber \\
&\stackrel{\text{(B.5b),(B.6a)}}{\le} 1 + \mathcal{C}_M^i \|u_{\ell}^{i+1} - u_{\ell}^i\|_{L^{\infty}(\Omega)}. \tag{B.8}
\end{align}
$$

The proof of the lower estimate is similar and uses that

$$
\|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}^i} \le (1 + \mathcal{C}_m^i \|u_{\ell}^{i+1} - u_{\ell}^i\|_{L^{\infty}(\Omega)})^{\frac{1}{2}} \|\mathcal{R}(\tilde{u})\|_{-1, u_{\ell}}.
$$

Finally, to see that the constants $\mathcal{C}_m^i$ are uniformly bounded, note that $\|u_\ell^i\|_{L^\infty(\Omega)}$ are bounded uniformly due to (B.6a), and consequently, $I_\ell^i(\mathbf{x})$ are subsets of a compact set $I_\ell \subset \mathbb{R}$ which does not depend on the iteration index $i$. The positivity and continuity of $F$ then implies that $\inf_{v \in I_\ell} F(v) =: F_m > 0$. Hence, we get

$$
\mathcal{C}_m^i \stackrel{(\text{B.5a})}{\le} \frac{1}{1-\alpha} \sup_{\mathbf{x} \in \Omega} \left\{ \sup_{v \in I_\ell^i(\mathbf{x})} \left| \frac{F(\mathbf{x},v) - F(\mathbf{x},u_\ell^i)}{v-u_\ell^i} \right| \frac{1}{\inf_{v \in I_\ell^i(\mathbf{x})} F(\mathbf{x},v)} \right\} \stackrel{(\text{B.1})}{\le} \frac{F_M}{(1-\alpha)F_m}. \quad (\text{B.9a})
$$

The same argument shows that also

$$
\mathcal{C}_M^i \le \frac{F_M}{(1 - \alpha)F_M}. \qquad (B.9b)
$$

□

The following result shows that if the iterations $\{u_{\ell}^i\}_{i \ge 0}$ are converging, then the iteration-dependent dual norm $\|\mathcal{R}(u_{\ell}^i)\|_{-1, u_{\ell}^i}$ converges to the error of the finite element solution measured in the iteration-independent natural discrete dual norm $\|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}}$.

**Corollary B.6** (Convergence of $\|\mathcal{R}(u_{\ell}^i)\|_{-1, u_{\ell}^i}$ to $\|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}}$). Let $\mathcal{R}: H_0^1(\Omega) \to H^{-1}(\Omega)$ from (1.1) be a continuous operator, and let the assumptions of Theorem B.5 hold. Then

$$
\|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} \to \|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}} \text{ as } i \to \infty.
$$

*Proof.* Using the triangle inequality,

$$
\left| \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} - \|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}} \right| \le \left| \|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}^{i}} - \|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}} \right| + \left| \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} - \|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}^{i}} \right|.
$$

The first term vanishes due to Theorem B.5, in particular the uniform bounds on $\mathcal{C}_m^i$ and $\mathcal{C}_M^i$ from (B.9). The second term on the right-hand-side vanishes due to the continuity of the operator $\mathcal{R}$. More precisely, as $u_{\ell}^i \to u_{\ell}$,

$$
\begin{align*}
\left| \|\mathcal{R}(u_{\ell}^i)\|_{-1, u_{\ell}^i} - \|\mathcal{R}(u_{\ell})\|_{-1, u_{\ell}^i} \right| &\le \|\mathcal{R}(u_{\ell}^i) - \mathcal{R}(u_{\ell})\|_{-1, u_{\ell}^i} \\
&= \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(u_{\ell}^i) - \mathcal{R}(u_{\ell}), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}^i}} \\
&\stackrel{\text{(B.2)}}{\le} \max \left( \sup_{\Omega} \sqrt{\frac{F(u_{\ell})}{F(u_{\ell}^i)}}, 1 \right) \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \mathcal{R}(u_{\ell}^i) - \mathcal{R}(u_{\ell}), \varphi \rangle}{\|\varphi\|_{1, u_{\ell}}} \\
&\le C \|\mathcal{R}(u_{\ell}^i) - \mathcal{R}(u_{\ell})\|_{-1, u_{\ell}} \to 0.
\end{align*}
$$

Here, the uniform boundedness of $F(u_\ell)/F(u_\ell^i)$ in $\Omega$ with respect to the iteration index $i$ has been used which can be proven following the arguments from Theorem B.5, see (B.8)–(B.9). $\square$

**Remark B.7** (Convergence in the iteration-dependent norm $\|\cdot\|_{1, u_{\ell}^i}$). Following Lemma 5.2, the statement of Corollary B.6 can be rewritten as

$$
\|u_{\ell}^{i} - u_{\langle \ell \rangle}^{i+1}\|_{1, u_{\ell}^{i}} \to \|u_{\ell} - u_{\langle \ell \rangle}\|_{1, u_{\ell}} \text{ as } i \to \infty.
$$
# Appendix C Error estimates for the Newton linearization for the gradient-independent diffusivity

In this section, we consider the Newton scheme for the gradient-independent diffusivity case (3.3). This leads to a nonsymmetric linearization problem which does not fit the symmetric framework of Sections 4 and 5.

## C.1 Newton iterations

The Newton(Raphson) scheme for iterative linearization of (4.3) in the gradient-independent diffusivity case (3.3), see [22], is defined in the spirit of (4.4), see also (A.2), but with an advection-reaction-diffusion linearization bilinear form $\mathcal{L}_N$: for $u_\ell^i \in V_\ell$ and $\xi, \zeta \in H_0^1(\Omega)$,

$$ \mathcal{L}_N(u_\ell^i; \xi, \zeta) := \langle L_N^i \xi, \zeta \rangle + (a_N^i \nabla \xi + w_N^i \xi, \nabla \zeta), \quad (C.1) $$

where $L_N^i$ and $a_N^i$ satisfy Assumption 4.1 taken for $L^i$ and $a^i$. In particular, with reference to (3.3) and Table 2,

$$ \begin{aligned} L_N^i(x) &:= \partial_\xi f(x, u_\ell^i), \quad a_N^i(x) := \tau \bar{K}(x) \mathcal{D}(x, u_\ell^i), \\ w_N^i(x) &:= \tau \bar{K}(x) [\partial_\xi \mathcal{D}(x, u_\ell^i) \nabla u_\ell^i + \partial_\xi q(x, u_\ell^i)]. \end{aligned} \quad (C.2) $$

Starting from an initial guess $u_\ell^0 \in V_\ell$, the Newton iterates $\{u_\ell^i\}_{i \ge 1} \subset V_\ell$ are then obtained by solving, similar to (4.7), the problems

$$ \mathcal{L}_N(u_\ell^i; u_\ell^{i+1} - u_\ell^i, \varphi_\ell) = -\langle \mathcal{R}(u_\ell^i), \varphi_\ell \rangle \quad \forall \varphi_\ell \in V_\ell. \quad (C.3) $$

Although $\mathcal{L}_N$ is bilinear, it is not symmetric, and thus it does not correspond to an inner product, in contrast to (4.4). Consequently, the orthogonality relation obtained in Theorem 5.3 does not hold in this case. However, using the techniques of Sections 5 and 8, it is still possible to provide an upper and a lower bound on an iteration-dependent dual norm of $\mathcal{R}(u_\ell^i)$ which converge to an orthogonal relation when the asymmetry (advection) vanishes. The following restrictions on $w_N^i$ are assumed for this purpose.

**Assumption C.1** (Properties of the advection vector $w_N^i$). For $\{u_\ell^i\}_{i \ge 0} \subset V_\ell$ stemming from the iterative Newton linearization (C.3), we assume that $w_N^i \in L^\infty(\Omega; \mathbb{R}^d)$ and there exists a constant $C_N \in [0, 2)$ such that

$$ w_N^i(x) (a_N^i(x))^{-1} w_N^i(x) \le C_N^2 L_N^i(x) \quad \text{for a.e. } x \in \Omega, \forall i \ge 0. \quad (C.4) $$

**Remark C.2** (Validity of Assumption C.1). For problems having linear diffusion coefficients, i.e. $\partial_\xi \mathcal{D} = 0$ in (3.3), Assumption C.1 holds if $L_N^i > 0$ and $\tau$ (which corresponds to the time-step size for time-discretized unsteady problems, cf. Class II in Example 3.4) is small enough. For problems having nonlinear diffusion coefficients, such as the Richards equation (see Section 6.2), the condition is still satisfied if $u_\ell^i$ is bounded in $W^{1,\infty}(\Omega)$ uniformly with respect to $i$ and the time-step size $\tau > 0$ is small. The time-step and non-degeneracy restrictions are well-known constraints in proving the convergence of the Newton scheme in applications, see [52] and the references therein.

Similarly to (4.15) and (4.18) we define, for any $i \ge 0$, the energy norms

$$ \|\xi\|_{1, u_\ell^i} := (\|(L_N^i)^{\frac{1}{2}} \xi\|^2 + \|(a_N^i)^{\frac{1}{2}} \nabla \xi\|^2)^{\frac{1}{2}}, \quad \|\zeta\|_{-1, u_\ell^i} := \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \zeta, \varphi \rangle}{\|\varphi\|_{1, u_\ell^i}}. \quad (C.5) $$

Contrary to (4.15), however, we do not have $\|\xi\|_{1, u_\ell^i}^2 = \mathcal{L}_N(u_\ell^i; \xi, \xi)$ for $\xi \in H_0^1(\Omega)$. Instead, by the Assumption C.1, there holds
**Lemma C.3** (Inequalities for $\mathcal{L}_N$ and the $||\cdot||_{\pm 1, u_\ell^i}$ norms). For all $\xi, \zeta \in H_0^1(\Omega)$,

$$|\mathcal{L}_N(u_\ell^i; \xi, \zeta)| \le (1 + C_N) ||\xi||_{1, u_\ell^i} ||\zeta||_{1, u_\ell^i}, \qquad (C.6a)$$

$$(1 - \frac{C_N}{2}) ||\xi||_{1, u_\ell^i}^2 \le \mathcal{L}_N(u_\ell^i; \xi, \xi) \le (1 + \frac{C_N}{2}) ||\xi||_{1, u_\ell^i}^2. \qquad (C.6b)$$

*Proof.* Observe from (C.1) and (C.5) that

$$|\mathcal{L}_N(u_\ell^i; \xi, \zeta)| \le ||\xi||_{1, u_\ell^i} ||\zeta||_{1, u_\ell^i} + |(w_N^i \xi, \nabla \zeta)|.$$

For (C.6a), we estimate the last term using the Cauchy-Schwarz inequality:

$$|(w_N^i \xi, \nabla \zeta)| \le \|(a_N^i)^{-\frac{1}{2}} w_N^i \xi\| \|(a_N^i)^{\frac{1}{2}} \nabla \zeta\| \qquad (C.4) \\ \le C_N \|(L_N^i)^{\frac{1}{2}} \xi\| \|(a_N^i)^{\frac{1}{2}} \nabla \zeta\| \le C_N \|\xi\|_{1, u_\ell^i} \|\zeta\|_{1, u_\ell^i}. \qquad (C.5)$$

For proving (C.6b), we note that

$$|(w_N^i \xi, \nabla \xi)| \stackrel{(C.4)}{\le} C_N \|(L_N^i)^{\frac{1}{2}} \xi\| \|(a_N^i)^{\frac{1}{2}} \nabla \xi\| \stackrel{(2.2)}{\le} \frac{C_N}{2} \left( \|(L_N^i)^{\frac{1}{2}} \xi\|^2 + \|(a_N^i)^{\frac{1}{2}} \nabla \xi\|^2 \right) = \frac{C_N}{2} \|\xi\|_{1, u_\ell^i}^2. \quad \square$$

## C.2 Bounds on the dual norm of the residual by linearization and discretization components

Extending the definition of $\mathcal{R}_{\text{disc}}^{u_\ell^i} : H_0^1(\Omega) \to H^{-1}(\Omega)$ from (5.1) to the Newton case, one has

$$\langle \mathcal{R}_{\text{disc}}^{u_\ell^i}(\xi), \varphi \rangle := \mathcal{L}_N(u_\ell^i; \xi - u_\ell^i, \varphi) + \langle \mathcal{R}(u_\ell^i), \varphi \rangle, \quad \varphi \in H_0^1(\Omega), \qquad (C.8)$$

so that using (C.3), $\langle \mathcal{R}_{\text{disc}}^{u_\ell^i}(u_\ell^{i+1}), \varphi_\ell \rangle = 0$ for all $\varphi_\ell \in V_\ell$. Then, in extension of Theorem 5.3, we have:

**Theorem C.4** (Upper and lower bounds on error for non-symmetric $\mathcal{L}_N$ (Newton)). In the context of the gradient-independent diffusivity case (3.3), let the bilinear form $\mathcal{L}_N$ be defined by (C.1), the sequence $\{u_\ell^i\}_{i \ge 0} \subset V_\ell$ by (C.3), the norms $||\cdot||_{\pm 1, u_\ell^i}$ by (C.5), and $\mathcal{R}_{\text{disc}}^{u_\ell^i} : H_0^1(\Omega) \to H^{-1}(\Omega)$ by (C.8). Let Assumption C.1 hold for some $C_N \in [0, 2)$. Then, one has

$$(1 - \frac{C_N}{2})^2 \left[ \frac{1}{(1+C_N)^2} \left\| \mathcal{R}_{\text{disc}}^{u_\ell^i}(u_\ell^{i+1}) \right\|_{-1, u_\ell^i}^2 + \left\| u_\ell^i - u_\ell^{i+1} \right\|_{1, u_\ell^i}^2 \right. \\ \left. - \frac{4C_N}{2-C_N} \left\| \mathcal{R}_{\text{disc}}^{u_\ell^i}(u_\ell^{i+1}) \right\|_{-1, u_\ell^i}^2 \left\| u_\ell^i - u_\ell^{i+1} \right\|_{1, u_\ell^i} \right] \\ \le \left\| \mathcal{R}(u_\ell^i) \right\|_{-1, u_\ell^i}^2 \le (1 + C_N)^2 \left[ \frac{4}{(2-C_N)^2} \left\| \mathcal{R}_{\text{disc}}^{u_\ell^i}(u_\ell^{i+1}) \right\|_{-1, u_\ell^i}^2 + \left\| u_\ell^i - u_\ell^{i+1} \right\|_{1, u_\ell^i}^2 \right. \\ \left. + \frac{4C_N}{2-C_N} \left\| \mathcal{R}_{\text{disc}}^{u_\ell^i}(u_\ell^{i+1}) \right\|_{-1, u_\ell^i}^2 \left\| u_\ell^i - u_\ell^{i+1} \right\|_{1, u_\ell^i} \right].$$

**Remark C.5** (Relation to Theorem 5.3). If the asymmetry (the advection term $w_N^i$) in $\mathcal{L}_N$ disappears, then the constant $C_N$ in (C.4) can be taken as 0. In this case, both the upper and the lower estimates of $||\mathcal{R}(u_\ell^i)||_{-1, u_\ell^i}$ in Theorem C.4 change to the orthogonal relation of Theorem 5.3.
Proof. Similarly to (4.13), let $u_{\langle\ell\rangle}^{i+1} \in H_0^1(\Omega)$ solve

$$ \mathcal{L}_N(u_{\langle\ell\rangle}^i; u_{\langle\ell\rangle}^{i+1} - u_{\ell}^i, \varphi) = -\langle \mathcal{R}(u_{\ell}^i), \varphi \rangle \quad \forall \varphi \in H_0^1(\Omega). \quad (C.9) $$

Then, following (C.3)–(C.9) and Lemma C.3, we have that

$$ \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}} \stackrel{(C.9)}{=} \sup_{\varphi \in H_{0}^{1}(\Omega)} \frac{\mathcal{L}_{N}(u_{\ell}^{i}; u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i}, \varphi)}{\|\varphi\|_{1, u_{\ell}^{i}}} \begin{cases} \stackrel{(C.6a)}{\le} (1 + C_{N}) \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i}\|_{1, u_{\ell}^{i}}, \\ \stackrel{(C.6b)}{\ge} (1 - \frac{C_{N}}{2}) \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i}\|_{1, u_{\ell}^{i}}. \end{cases} \quad (C.10) $$

Expanding $\|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i}\|_{1, u_{\ell}^{i}}^2$, one has

$$ \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i}\|_{1, u_{\ell}^{i}}^{2} \stackrel{(C.5)}{=} \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}}^{2} + \|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{1, u_{\ell}^{i}}^{2} \\ + 2[(L_{N}^{i}(u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}), u_{\ell}^{i+1} - u_{\ell}^{i}) + (a_{N}^{i}\nabla(u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}), \nabla(u_{\ell}^{i+1} - u_{\ell}^{i}))]. \quad (C.11a) $$

Observe that since $u_{\ell}^{i+1} - u_{\ell}^{i} \in V_{\ell}$, we have by the definition of finite element solutions (C.3) that

$$ \mathcal{L}_N(u_{\ell}^{i}; u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}, u_{\ell}^{i+1} - u_{\ell}^{i}) = 0. \quad (C.11b) $$

Hence, we have following the proof of Lemma C.3,

$$ \begin{aligned} & |(L_N^i(u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}), u_{\ell}^{i+1} - u_{\ell}^{i}) + (a_N^i \nabla(u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}), \nabla(u_{\ell}^{i+1} - u_{\ell}^{i}))| \\ & \stackrel{(C.1)}{=} |\mathcal{L}_N(u_{\ell}^{i}; u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}, u_{\ell}^{i+1} - u_{\ell}^{i}) - (w_N^i(u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}), \nabla(u_{\ell}^{i+1} - u_{\ell}^{i}))| \quad (C.11c) \\ & \stackrel{(C.11b)}{=} |(w_N^i(u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}), \nabla(u_{\ell}^{i+1} - u_{\ell}^{i}))| \stackrel{(C.7)}{\le} C_N \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}} \|u_{\ell}^{i+1} - u_{\ell}^{i}\|_{1, u_{\ell}^{i}}. \end{aligned} $$

Finally, noting that

$$ \begin{aligned} \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}} &\stackrel{(C.6b)}{\le} \frac{2}{2 - C_N} \frac{\mathcal{L}_N(u_{\ell}^{i}; u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}, u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1})}{\|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}}} \\ &\stackrel{(C.8),(C.9)}{\le} \sup_{\varphi \in H_0^1(\Omega)} \frac{\mathcal{L}_N(u_{\ell}^{i}; u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}, \varphi)}{\|\varphi\|_{1, u_{\ell}^{i}}} \frac{2}{2 - C_N} \|\mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})\|_{-1, u_{\ell}^{i}}, \end{aligned} \quad (C.11d) $$

we have the upper bound. The proof of the lower bound is similar, but requires the estimation

$$ \|\mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})\|_{-1, u_{\ell}^{i}} \le (1 + C_N) \|u_{\langle\ell\rangle}^{i+1} - u_{\ell}^{i+1}\|_{1, u_{\ell}^{i}}, $$

which follows from the last equality in (C.11d) together with (C.6a). $\square$

### C.3 A posteriori error bounds

With the above developments, the a posteriori analysis of Sections 5 and 8 carries over to the Newton scheme. In particular, the counterpart of Theorem 5.7 for the Newton scheme becomes:

**Theorem C.6** (Guaranteed, efficient, and robust a posteriori estimates for non-symmetric $\mathcal{L}_N$ (Newton)). Let the assumptions of Theorem C.4 be satisfied. Rearranging (C.3), write it in the form (4.9) with $S^i$ given by (4.12a) and $\mathcal{F}^i := \tau \bar{\mathbf{K}}(\cdot) \mathbf{q}(\cdot, u_{\ell}^i) + (\tau \bar{\mathbf{K}}(\cdot) \mathcal{D}(\cdot, u_{\ell}^i) - \mathbf{a}_N^i) \nabla u_{\ell}^i +$
$w_N^i(u_\ell^{i+1} - u_\ell^i) \stackrel{(C.2)}{=} \tau \bar{K}(\cdot) q(\cdot, u_\ell^i) + w_N^i(u_\ell^{i+1} - u_\ell^i)$, which only has the reaction-diffusion operator given by $L_N^i$ and $a_N^i$ on the left-hand side and the computed $u_\ell^{i+1}$ on the right-hand side. Define a posteriori estimators as in (5.12), following the steps of Section 5.2 with in particular the equilibration given by Definition 5.4. Let $\vartheta_{\Omega, \ell, i} \ge 1$ be as in Theorem 5.7. Finally, let the constants $\theta_{N,m}$ and $\theta_{N,M}$, which converge to 1 as $C_N \to 0$, be defined as

$$ \theta_{N,m} := \inf_{t \ge 0} \frac{t^2 - C_N(2-C_N)t + \frac{1}{(1+C_N)^2}}{(1-C_N/2)^{-2}(t^2+1)} > -\infty, \quad \theta_{N,m} := \sup_{t \ge 0} \frac{t^2 + \frac{4C_N}{2-C_N}t + \frac{4}{(2-C_N)^2}}{(1+C_N)^{-2}(t^2+1)} \in (0, \infty). \tag{C.12} $$

Then, there holds

$$ \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2} \le \theta_{N,M} [\eta_{\Omega}^{i}]^{2}, \qquad (\text{guaranteed reliability}) $$

and, assuming $\theta_{N,m} > 0$,

$$ \theta_{N,m}[\eta_{\Omega}^{i}]^{2} \lesssim \vartheta_{\Omega, \ell, i}^{2} \|\mathcal{R}(u_{\ell}^{i})\|_{-1, u_{\ell}^{i}}^{2} + \theta_{N,m} [[\eta_{\text{osc},\Omega}^{i}]^{2} + \vartheta_{\Omega, \ell, i}^{2}[\eta_{\text{quad,F},\Omega}^{i}]^{2} + [\eta_{\text{quad,S},\Omega}^{i}]^{2} + [\eta_{\text{quad,S,osc},\Omega}^{i}]^{2}], \tag{global efficiency} $$

where the hidden constant in $\lesssim$ has the same dependence as in Theorem 5.7.

**Remark C.7** (The asymmetry factors from (C.12) $\theta_{N,m}$ and $\theta_{N,M}$). The asymmetry factors $\theta_{N,m}$ and $\theta_{N,M}$ are introduced to simplify in Theorem C.6 the bounds introduced in Theorem C.4 and to make them similar to those of Theorem 5.7. Observe that $\theta_{N,m}$ and $\theta_{N,M}$ are computable by calculating the extremas of the functions in (C.12). Thus, similarly to the variability constant $\vartheta_{\Omega, \ell, i}$, they can be estimated a posteriori.

*Proof.* As in Section 8, the dual norm of the discretization residual $\mathcal{R}_{\text{disc}}^{u_{\ell}^{i}}(u_{\ell}^{i+1})$ is estimated above and below by (8.10) and (5.14b) respectively, since the steps to obtain the estimates remain exactly the same after redefining $\mathcal{F}^i$. Then the estimates are proven by using (C.12) in Theorem C.4. We omit the details for the sake of conciseness. $\square$

## C.4 Numerical experiments

We repeat here the numerical experiments of Section 6.2 for the Newton linearization (C.3). The effectivity indices are defined as in (6.2), where $u_{\text{ref}}^{i+1} \in V_{10}$ now solves the problem

$$ (L_N^i(u_{\text{ref}}^{i+1} - u_\ell^i), \varphi_{\text{ref}}) + (a_N^i \nabla(u_{\text{ref}}^{i+1} - u_\ell^i), \nabla \varphi_{\text{ref}}) = -\langle \mathcal{R}(u_\ell^i), \varphi_{\text{ref}} \rangle \quad \forall \varphi_{\text{ref}} \in V_{10}, \tag{C.13} $$

thus, yielding $\|u_\ell^i - u_{\text{ref}}^{i+1}\|_{1, u_\ell^i} \approx \|\mathcal{R}(u_\ell^i)\|_{-1, u_\ell^i}$. Figure 14 shows the results. The global effectivity indices defined in (6.2) can go below 1 in this case, since the factor $\theta_{N,M}$ is not included. This is observed for $\tau = 1$ at iteration $i = 3$. However, the global effectivity indices stay between 1 and 2.5 in all other cases. Local effectivity indices also remain close to 1 in most of the domain. We find the agreement to be excellent. However, it is to be noted that the estimators only indicate the actual error if Assumption C.1 holds, which cannot always be guaranteed a priori. The quadratic convergence property of the Newton scheme is exhibited in the right plot of Figure 14.
<figure><img src="image_16.png" /><figcaption>Figure 14: [Appendix C, setting of Section 6.2, $\tau$ varies] Numerical results for the Newton scheme, $\tau = 1$ (top) and $\tau = 0.01$ (bottom). Global effectivity indices (left); local effectivity indices for $\ell = 2$ at the iteration $i = \bar{i}$ from (5.16) (center); the total error $\|\mathcal{R}(u_{\ell}^i)\|_{-1, u_{\ell}^i}$, the linearization estimator $\eta_{\text{lin},\Omega}^i$, and the total estimator $\eta_{\Omega}^i$ against the iteration index $i$ (right, logarithmic scale). The vertical black-dotted lines indicate the stopping criterion $i = \bar{i}$ from (5.16) with $\mu = 0.05$. The quadratic convergence property of the Newton scheme can clearly be observed from the ever-steepening slope of $\eta_{\text{lin},\Omega}^i$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: [Appendix C, setting of Section 6.2, $\tau$ varies] Global effectivity indices (6.2a) vs. $1/\tau$ for the Newton linearization at different discretization levels $\ell$ for $i = \bar{i}$ from (5.16). The variability constant $\vartheta_{\Omega,\ell,i}$ from (5.15) is also plotted (scaled by half and for $\ell = 2$).</figcaption></figure>

Figure 15 plots the global effectivity indices (6.2a) against the strength of the nonlinearity represented by $1/\tau$. These indices remain between 1.5 and 2.6, numerically confirming the robustness of our a posteriori estimators. In fact, the obtained effectivity indices for the Newton scheme are almost identical to those of the Picard scheme shown in Figure 9. This results from the fact that the advection term $(w_N^i(u_\ell^{i+1} - u_\ell^i), \nabla\varphi_\ell)$ in (C.3) is negligible at $i = \bar{i}$ compared to the symmetric part of $\mathcal{L}_N$, which makes the two schemes almost coincide.
# Acknowledgements

K. Mitra was supported by Fonds voor Wetenschappelijk Onderzoek (FWO) through the Junior Postdoctoral Fellowship during the completion of this work.

# References

[1] Ainsworth, M., and Babuška, I. Reliable and robust a posteriori error estimation for singularly perturbed reaction-diffusion problems. *SIAM J. Numer. Anal.* **36** (1999), 331–353.

[2] Ainsworth, M., and Vejchodský, T. Robust error bounds for finite element approximation of reaction-diffusion problems with non-constant reaction coefficient in arbitrary space dimension. *Comput. Methods Appl. Mech. Engrg.* **281** (2014), 184–199. http://dx.doi.org/10.1016/j.cma.2014.08.005.

[3] Bartels, S., and Milicevic, M. Primal-dual gap estimators for a *posteriori* error analysis of nonsmooth minimization problems. *ESAIM Math. Model. Numer. Anal.* **54** (2020), 1635–1660. https://doi.org/10.1051/m2an/2019074.

[4] Bernardi, C., Dakroub, J., Mansour, G., and Sayah, T. A posteriori analysis of iterative algorithms for a nonlinear problem. *J. Sci. Comput.* **65** (2015), 672–697. https://doi.org/10.1007/s10915-014-9980-4.

[5] Bernardi, C., and Verfürth, R. Adaptive finite element methods for elliptic equations with non-smooth coefficients. *Numer. Math.* **85** (2000), 579–608. https://doi.org/10.1007/PL00005393.

[6] Blechta, J., Málek, J., and Vohralík, M. Localization of the $W^{-1,q}$ norm for local a posteriori efficiency. *IMA J. Numer. Anal.* **40** (2020), 914–950. https://doi.org/10.1093/imanum/drz002.

[7] Boffi, D., Brezzi, F., and Fortin, M. *Mixed finite element methods and applications*, vol. **44** of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013. https://doi.org/10.1007/978-3-642-36519-5.

[8] Botti, M., and Riedlbeck, R. Equilibrated stress tensor reconstruction and a posteriori error estimation for nonlinear elasticity. *Comput. Methods Appl. Math.* **20** (2020), 39–59. https://doi.org/10.1515/cmam-2018-0012.

[9] Braess, D., Pillwein, V., and Schöberl, J. Equilibrated residual error estimates are p-robust. *Comput. Methods Appl. Mech. Engrg.* **198** (2009), 1189–1197. http://dx.doi.org/10.1016/j.cma.2008.12.010.

[10] Braess, D., and Schöberl, J. Equilibrated residual error estimator for edge elements. *Math. Comp.* **77** (2008), 651–672. http://dx.doi.org/10.1090/S0025-5718-07-02080-7.

[11] Buhr, A., Engwer, C., Ohlberger, M., and Rave, S. ArbiLoMod, a simulation technique designed for arbitrary local modifications. *SIAM J. Sci. Comput.* **39** (2017), A1435–A1465. https://doi.org/10.1137/15M1054213.

[12] Cai, D., Cai, Z., and Zhang, S. Robust equilibrated a posteriori error estimator for higher order finite element approximations to diffusion problems. *Numer. Math.* **144** (2020), 1–21. https://doi.org/10.1007/s00211-019-01075-1.

[13] Celia, M., Bouloutas, E., and Zarba, R. General mass-conservative numerical solution for the unsaturated flow equation. *Water Resources Research* **26** (1990), 1483–1496.

[14] Chaillou, A., and Suri, M. A posteriori estimation of the linearization error for strongly monotone nonlinear operators. *J. Comput. Appl. Math.* **205** (2007), 72–87. http://dx.doi.org/10.1016/j.cam.2006.04.041.

[15] Chaillou, A. L., and Suri, M. Computable error estimators for the approximation of non-linear problems by linearized models. *Comput. Methods Appl. Mech. Engrg.* **196** (2006), 210–224. http://dx.doi.org/10.1016/j.cma.2006.03.008.
[16] Ciarlet, Jr., P., and Vohralík, M. Localization of global norms and robust a posteriori error control for transmission problems with sign-changing coefficients. *ESAIM Math. Model. Numer. Anal.* **52** (2018), 2037–2064. https://doi.org/10.1051/m2an/2018034.

[17] Ciarlet, P. G. *The Finite Element Method for Elliptic Problems*, vol. **4** of *Studies in Mathematics and its Applications*. North-Holland, Amsterdam, 1978.

[18] Cohen, A., Dahmen, W., and DeVore, R. Adaptive wavelet schemes for nonlinear variational problems. *SIAM J. Numer. Anal.* **41** (2003), 1785–1823. http://dx.doi.org/10.1137/S0036142902412269.

[19] Cottereau, R., Chamoin, L., and Díez, P. Strict error bounds for linear and nonlinear solid mechanics problems using a patch-based flux-free method. *Mechanics & Industry* **11** (2010), 249–254. https://doi.org/10.1051/meca/2010049.

[20] Cottereau, R., Díez, P., and Huerta, A. Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method. *Comput. Mech.* **44** (2009), 533–547. https://doi.org/10.1007/s00466-009-0388-1.

[21] Destuynder, P., and Métivet, B. Explicit error bounds in a conforming finite element method. *Math. Comp.* **68** (1999), 1379–1396. http://dx.doi.org/10.1090/S0025-5718-99-01093-5.

[22] Deuflhard, P. *Newton methods for nonlinear problems: affine invariance and adaptive algorithms*, vol. **35**. Springer Science & Business Media, Berlin, 2005.

[23] El Alaoui, L., Ern, A., and Vohralík, M. Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems. *Comput. Methods Appl. Mech. Engrg.* **200** (2011), 2782–2795. http://dx.doi.org/10.1016/j.cma.2010.03.024.

[24] Ern, A., and Vohralík, M. Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs. *SIAM J. Sci. Comput.* **35** (2013), A1761–A1791. http://dx.doi.org/10.1137/120896918.

[25] Ern, A., and Vohralík, M. Stable broken $H^1$ and $H(\text{div})$ polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions. *Math. Comp.* **89** (2020), 551–594. http://dx.doi.org/10.1090/mcom/3482.

[26] Evans, L. *Partial differential equations*, vol. **19**. American Mathematical Society, 2022.

[27] Gantner, G., Haberl, A., Praetorius, D., and Schimanko, S. Rate optimality of adaptive finite element methods with respect to overall computational costs. *Math. Comp.* **90** (2021), 2011–2040. https://doi.org/10.1090/mcom/3654.

[28] Garau, E. M., Morin, P., and Zuppa, C. Convergence of an adaptive Kačanov FEM for quasi-linear problems. *Appl. Numer. Math.* **61** (2011), 512–529. https://doi.org/10.1016/j.apnum.2010.12.001.

[29] Guo, M., Han, W., and Zhong, H. Legendre-Fenchel duality and a generalized constitutive relation error. arXiv preprint 1611.05589, https://arxiv.org/abs/1611.05589, 2016.

[30] Haberl, A., Praetorius, D., Schimanko, S., and Vohralík, M. Convergence and quasi-optimal cost of adaptive algorithms for nonlinear operators including iterative linearization and algebraic solver. *Numer. Math.* **147** (2021), 679–725. https://doi.org/10.1007/s00211-021-01176-w.

[31] Han, W. A posteriori error analysis for linearization of nonlinear elliptic problems and their discretizations. *Math. Methods Appl. Sci.* **17** (1994), 487–508. http://dx.doi.org/10.1002/mma.1670170702.

[32] Han, W. *A posteriori error analysis via duality theory*, vol. **8** of *Advances in Mechanics and Mathematics*. Springer-Verlag, New York, 2005.
[33] Harnist, A., Mitra, K., Rappaport, A., and Vohralík, M. Robust energy a posteriori estimates for nonlinear elliptic problems. HAL Preprint 04033438, https://hal.inria.fr/hal-04033438, 2023.

[34] Hecht, F. New development in FreeFem++. *J. Numer. Math.* **20** (2012), 251–265. https://doi.org/10.1515/jnum-2012-0013.

[35] Heid, P., Praetorius, D., and Wihler, T. P. Energy contraction and optimal convergence of adaptive iterative linearized finite element methods. *Comput. Methods Appl. Math.* **21** (2021), 407–422. https://doi.org/10.1515/cmam-2021-0025.

[36] Heid, P., and Wihler, T. Adaptive iterative linearization Galerkin methods for nonlinear problems. *Mathematics of Computation* **89** (2020), 2707–2734. https://doi.org/10.1090/mcom/3545.

[37] Holst, M., Szypowski, R., and Zhu, Y. Adaptive finite element methods with inexact solvers for the nonlinear Poisson-Boltzmann equation. In *Domain decomposition methods in science and engineering XX*, vol. **91** of Lect. Notes Comput. Sci. Eng. Springer, Heidelberg, 2013, pp. 167–174. https://doi.org/10.1007/978-3-642-35275-1.

[38] Houston, P., Süli, E., and Wihler, T. A posteriori error analysis of hp-version discontinuous Galerkin finite-element methods for second-order quasi-linear elliptic PDEs. *IMA Journal of Numerical Analysis* **28** (2008), 245–273. https://doi.org/10.1093/imanum/drm009.

[39] Jäger, W., and Kačur, J. Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes. *ESAIM: Mathematical Modelling and Numerical Analysis* **29** (1995), 605–627.

[40] Karátson, J., and Korotov, S. Sharp upper global a posteriori error estimates for nonlinear elliptic variational problems. *Appl. Math.* **54** (2009), 297–336. http://dx.doi.org/10.1007/s10492-009-0020-x.

[41] Kim, K.-Y. A posteriori error estimators for locally conservative methods of nonlinear elliptic problems. *Appl. Numer. Math.* **57** (2007), 1065–1080. http://dx.doi.org/10.1016/j.apnum.2006.09.010.

[42] Kopteva, N. Energy-norm a posteriori error estimates for singularly perturbed reaction-diffusion problems on anisotropic meshes. *Numer. Math.* **137** (2017), 607–642. https://doi.org/10.1007/s00211-017-0889-3.

[43] Kopteva, N. Fully computable a posteriori error estimator using anisotropic flux equilibration on anisotropic meshes. ArXiv Preprint 1704.04404, https://arxiv.org/abs/1704.04404, 2017.

[44] Kovács, B. On the numerical performance of a sharp a posteriori error estimator for some nonlinear elliptic problems. *Appl. Math.* **59** (2014), 489–508. https://doi.org/10.1007/s10492-014-0068-0.

[45] Ladevèze, P., and Moës, N. A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control. In *Advances in adaptive computational methods in mechanics (Cachan, 1997)*, vol. **47** of Stud. Appl. Mech. Elsevier Sci. B. V., Amsterdam, 1998, pp. 231–256. https://doi.org/10.1016/S0922-5382(98)80013-5.

[46] Ladevèze, P., and Pelle, J.-P. *Mastering calculations in linear and nonlinear mechanics*. Springer-Verlag, New York, 2005.

[47] Ladyzhenskaya, O., Uraltseva, N., and Ehrenpreis, L. *Linear and quasilinear elliptic equations*. Academic Press, 1968.

[48] Li, H., Wu, Z., Yin, J., and Zhao, J. *Nonlinear diffusion equations*. World Scientific, Singapore, 2001.

[49] Lions, P. On the existence of positive solutions of semilinear elliptic equations. *SIAM review* **24** (1982), 441–467.
[50] List, F., and Radu, F. A study on iterative methods for solving Richards' equation. *Computational Geosciences* **20** (2016), 341–353. https://doi.org/10.1007/s10596-016-9566-3.

[51] Meyers, N. G. An $L^p$-estimate for the gradient of solutions of second order elliptic divergence equations. *Ann. Scuola Norm. Sup. Pisa Cl. Sci.* **(3)** **17** (1963), 189–206.

[52] Mitra, K., and Pop, I. A modified L-scheme to solve nonlinear diffusion problems. *Computers & Mathematics with Applications* **77** (2019), 1722–1738. 7th International Conference on Advanced Computational Methods in Engineering (ACOMEN 2017), https://doi.org/10.1016/j.camwa.2018.09.042.

[53] Mitra, K., and Vohralík, M. A posteriori error estimates for the Richards equation. HAL Preprint 03328944, submitted for publication, https://hal.inria.fr/hal-03328944, 2022.

[54] Neittaanmaki, P., and Repin, S. A posteriori error identities for nonlinear variational problems. *Ann. Acad. Rom. Sci. Ser. Math. Appl.* **7** (2015), 157–172.

[55] Pop, I., Radu, F., and Knabner, P. Mixed finite elements for the Richards equation: linearization procedure. *Journal of Computational and Applied Mathematics* **168** (2004), 365–373. https://doi.org/10.1016/j.cam.2003.04.008.

[56] Pousin, J., and Rappaz, J. Consistency, stability, a priori and a posteriori errors for Petrov-Galerkin methods applied to nonlinear problems. *Numer. Math.* **69** (1994), 213–231. http://dx.doi.org/10.1007/s002110050088.

[57] Repin, S. *A posteriori estimates for partial differential equations*, vol. **4** of *Radon Series on Computational and Applied Mathematics*. Walter de Gruyter GmbH & Co. KG, Berlin, 2008. http://dx.doi.org/10.1515/9783110203042.

[58] Repin, S. I. A posteriori error estimation for nonlinear variational problems by duality theory. *Zap. Nauchn. Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov. (POMI)* **243** (1997), 201–214, 342.

[59] Repin, S. I. A posteriori error estimation for variational problems with uniformly convex functionals. *Math. Comp.* **69** (2000), 481–500. http://dx.doi.org/10.1090/S0025-5718-99-01190-4.

[60] Smears, I., and Vohralík, M. Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems. *ESAIM Math. Model. Numer. Anal.* **54** (2020), 1951–1973. https://doi.org/10.1051/m2an/2020034.

[61] Stokke, J. S., Mitra, K., Storvik, E., Both, J. W., and Radu, F. A. An adaptive solution strategy for Richards' equation. arXiv 2301.02055, https://doi.org/10.48550/arXiv.2301.02055, 2023.

[62] van Genuchten, M. A closed form for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Amer. J.* **44** (1980), 892–898.

[63] Verfürth, R. A posteriori error estimates for nonlinear problems. Finite element discretizations of elliptic equations. *Math. Comp.* **62** (1994), 445–475. http://www.jstor.org/stable/2153518.

[64] Verfürth, R. Robust a posteriori error estimators for a singularly perturbed reaction-diffusion equation. *Numer. Math.* **78** (1998), 479–493. https://doi.org/10.1007/s002110050322.

[65] Vohralík, M. Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients. *J. Sci. Comput.* **46** (2011), 397–438. http://dx.doi.org/10.1007/s10915-010-9410-1.

[66] Zarantonello, E. *Solving functional equations by contractive averaging*. Mathematics Research Center, United States Army, University of Wisconsin, 1960.
[67] Zeidler, E. *Applied Functional Analysis: Applications to Mathematical Physics (Applied Mathematical Sciences)* (v. 108). Springer, 1995.

[68] Zeidler, E. *Nonlinear Functional Analysis and Its Applications: II/B: Nonlinear Monotone Operators*. Springer Science & Business Media, 2013.