Finite element error estimates in non-energy norms for the two-dimensional scalar Signorini problem

This paper is concerned with error estimates for the piecewise linear finite element approximation of the two-dimensional scalar Signorini problem on a convex polygonal domain Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document}. Using a Céa-type lemma, a supercloseness result, and a non-standard duality argument, we prove W1,p(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W^{1,p}(\varOmega )$$\end{document}-, L∞(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^\infty (\varOmega )$$\end{document}-, W1,∞(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W^{1,\infty }(\varOmega )$$\end{document}-, and H1/2(∂Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^{1/2}(\partial \varOmega )$$\end{document}-error estimates under reasonable assumptions on the regularity of the exact solution and Lp(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^p(\varOmega )$$\end{document}-error estimates under comparatively mild assumptions on the involved contact sets. The obtained orders of convergence turn out to be optimal for problems with essentially bounded right-hand sides. Our results are accompanied by numerical experiments which confirm the theoretical findings.


Introduction
The aim of this paper is to study the accuracy of the piecewise linear finite element method for the two-dimensional scalar Signorini problem on a convex polygonal domain Ω (in our analysis, for simplicity, assumed to be the unit square-see Sect. 2 for the precise assumptions). Using a Céa-type lemma, a supercloseness result, and a non-standard duality argument that is based on ideas of B Constantin Christof christof@ma.tum.de 1 Mosco, we establish new W 1, p -and L p -error estimates for the problem (1) that, in view of the W 2, p -and H s -regularity properties of the exact solution u, are optimal for right-hand sides f ∈ L ∞ (Ω). In particular, we prove an L 4 -error estimate of the form u − u h L 4 (Ω) = O(h 2−ε ) for all ε > 0 which explains the order of convergence in the lower L p -norms that is typically observed in numerical experiments, cf. Sect. 7 and [41,Section 7]. For the main contributions of this paper, we refer the reader to the regularity result in Theorem 2.3, the supercloseness result in Theorem 3.6, and the a priori finite element error estimates collected in Theorems 4.3 and 6.1.
Before we begin with our analysis, let us give some background: As one of the simplest examples of a problem that models contact, the Signorini problem (1) (along with its various reformulations and the closely related obstacle problem) has been subject to active research for a long period of time. In the context of finite element methods, early contributions on (1) and its approximation can be traced back at least to the nineteen-seventies, see [8,9,24,33,37], and even though these seminal works have been followed by a large number of other papers, the analysis of, for instance, FE-error estimates for (1) still receives considerable attention to this day. See, e.g., [16,21,23,41] for some recent contributions. The main reason for the ongoing interest in the problem (1) and the fact that, even after more than forty years, the approximation properties of the finite element method for (1) are still not fully understood is that the weak formulation of (1) takes the form of an elliptic variational inequality. This causes the differences u − u h between the continuous solution u of (1) and its finite element approximations u h to lack the property of Galerkin orthogonality and renders standard tools for the error analysis of finite element methods inapplicable. We remark that a notable exception to this rule are L ∞ -estimates which can be established along roughly the same lines as for elliptic PDEs in the case of the Signorini problem (1) by employing the discrete maximum principle or regularization techniques, cf. [4,11,13,15,18,22,28,34,35].
For error estimates in the H 1 -norm, which, as the energy norm of the problem, is the most natural choice for the analysis of (1), the lack of Galerkin orthogonality proved a challenge that could not be properly overcome for a considerable amount of time. Compare, for instance, with the contributions [5,8,9] in this context, which all require additional assumptions on the structure of the contact set {x ∈ ∂Ω | u(x) = 0} of the exact solution u to establish H 1 -error estimates of optimal order for the Signorini problem. Only very recently, it has been shown in [16] that these conditions are, in fact, not needed and that the H 1 -norm of u − u h is indeed always of order O(h) when (1) is discretized with standard piecewise linear finite elements (i.e., as in Sect. 2.3) and u possesses the natural regularity u ∈ H 2 (Ω) (cf. Proposition 2.1). We remark that, for other discretization schemes as, e.g., primal-dual approaches or Nitsche's method, similar results on the H 1 -error have recently also been obtained in [21]. See [10] for a survey article on this topic.
While H 1 -and L ∞ -error estimates for obstacle-and Signorini-type problems have been discussed quite extensively by various authors, results on the finite element error in other norms are only very rarely addressed in the literature. See, e.g., [12,31,33,40,41,43] for some of the few contributions on this topic. The reason for this is that, for error estimates in general L p -and W 1, p -norms, the missing Galerkin orthogonality of the differences u − u h is an even more severe problem than in the H 1 -case. This becomes apparent, for instance, in the study of the L 2 -error: Recall that, for linear elliptic partial differential equations, L 2 -error estimates of optimal order are typically established by means of the so-called Aubin-Nitsche trick. This trick is based on the idea to consider a dual partial differential equation, which contains the primal error u − u h as a right-hand side, and requires three main ingredients: an error estimate of optimal order in the H 1 -norm, the Galerkin orthogonality of the finite element approximations, and the H 2 -regularity of the dual solution, see [6,Section 7], [14,Section 3.2] or other standard references. Extending such a duality argument to the Signorini problem (1), where Galerkin orthogonality is not available, is clearly far from trivial. Nevertheless, since the derivation of the first H 1 -error estimates for elliptic variational inequalities with unilateral constraints in the nineteen-seventies, several authors have tried to accomplish precisely that, see [31,33,40,43]. The approaches that have been proposed in this context are typically based on the idea to consider a suitably defined dual variational inequality that, due to its construction, allows to bypass the lack of Galerkin orthogonality of the primal error u − u h . Unfortunately, at least to the best of the authors' knowledge, none of the contributions published so far has been able to simultaneously also satisfy the third prerequisite of the Aubin-Nitsche trick, namely, to show that the dual solution possesses enough regularity for the classical duality argument to go through. Compare, e.g., with [33] and [43,Section 5.2] in this context, where the H 2 -regularity of the dual solution is used as an assumption, or with [40] where it is implicitly assumed that the dual problem possesses a sufficiently regular solution, satisfies a constraint qualification, and admits sufficiently regular multipliers. Approaches that follow different lines to establish L 2 -error estimates, on the other hand, typically yield orders of convergence that are far from optimal. See, e.g., [41], where an error estimate in the H 1/2 -norm is derived on the boundary and subsequently used to establish an L 2 -estimate of order O(h 3/2−ε ) for all ε > 0. The results on the L 2 -error (and the error in the lower L p -norms in general) available in the literature are thus not very satisfactory and prove to be very unnatural in view of the required regularity assumptions on the involved primal and dual quantities. Even less seems to be known about W 1, p -error estimates for the Signorini problem with p = 2. At least to the authors' best knowledge, there are no contributions on this topic.
The purpose of the present paper is to demonstrate that, in the situation of (1), it is indeed possible to derive finite element error estimates of optimal order in non-energy norms for the Signorini problem while working only with reasonable assumptions on primal quantities. To be more precise, in what follows, we show that, if u enjoys a composite W 2, p -and H s -regularity, that can be proved to hold in various situations (see Theorem 2.3), then it is possible to establish error estimates in W 1, p (Ω), L ∞ (Ω), W 1,∞ (Ω), and H 1/2 (∂Ω) for (1) that are optimal for problems with L ∞ (Ω)right-hand sides (see Corollaries 4.1 and 4.2 and Theorem 4.3). Under the additional assumption that the contact set of the continuous solution u is sufficiently regular and that the contact sets of the finite element approximations u h do not exhibit a degenerative behavior in the limit h 0 (and thus in a setting comparable to those in [5,8,9] -see conditions (A) and (A h ) in Sects. 2 and 5), we are moreover able to extend the classical Aubin-Nitsche trick to (1) and to establish an L 4 -estimate of the form u − u h L 4 (Ω) = O(h 2−ε ) for all ε > 0. In combination with our previous findings, this provides us with a complete set of optimal-order finite element error estimates for problems (1) with L ∞ (Ω)-right-hand sides that does not require any artificial assumptions on the regularity properties of dual quantities (see Theorems 4.3 and 6.1).
The method of proof that we use for our finite element error analysis is somewhat non-standard in that it does not rely on a multiplier reformulation of (1) but on certain one-sided approximation results that are apparently only rarely employed in the literature. For a problem on the unit square Ω = (0, 1) 2 , whose right-hand side f is in L ∞ (Ω) and whose solution u has a sufficiently regular contact set, our approach can essentially be summarized as follows: Using an elementary argument, we prove that the H 1 -error between the finite element approximation u h and the Ritz projection R h (u) of the exact solution u is smaller than the H 1 -norm of every finite element function w h that satisfies Lemma 3.4). This best approximation property yields, in tandem with results on unilateral finite element approximations (see Lemma 3.5) and the regularity properties of solutions to (1) Theorem 3.6). By exploiting this supercloseness property, inverse estimates, and standard results for the Ritz projection, we immediately arrive at error estimates of optimal order in W 1,4 (Ω), W 1,∞ (Ω), L ∞ (Ω), and H 1/2 (∂Ω) (see Theorem 4.3 and Remark 4.4). To study the error in the lower L p -norms, we follow an approach of Mosco and consider two dual problems, one for each of the components max(0, u − u h ) and min(0, u − u h ) (see Sect. 5). As we will see, the solutions of our dual variational inequalities suffer from the same regularity problems as those in [31,33,43] and cannot be expected to be elements of H 2 (Ω). However, by invoking the results of [19,20], we can show that W 2,4/3−ε -regularity for all ε > 0 is obtainable instead. This observation and the fact that q := 4/3 is precisely the Hölder conjugate of p := 4 allow us to invoke our W 1,4 -estimate to compensate the lack of regularity of the dual solutions and to arrive at an estimate of the type u − u h L 4 (Ω) = O(h 2−ε ) for all ε > 0 (see Theorem 6.1).
We would like to point out that the H 1/2 (∂Ω)-error estimate that we establish in Theorem 4.3 reproduces [41, Theorem 2.2] under slightly different assumptions on the regularity of the exact solution u (or the right-hand side f , respectively). Further, Theorem 4.3 shows that the order of convergence 3/2−ε that has been obtained in [41,Corollary 5.8] in L 2 (Ω) is, in fact, even achieved in the L ∞ -norm. Surprisingly, we obtain this L ∞ -result without ever invoking the discrete maximum principle (which is normally used to prove pointwise error estimates for variational inequalities with unilateral constraints) and without the related assumptions on the underlying triangulation, cf. [4,11,13,15,18,28,34]. Theorem 6.1 finally improves the order of convergence in [41,Corollary 5.8] by the factor h 1/2 and yields an L 4 -error estimate that is optimal. To the best of our knowledge, the L p -and W 1, p -error estimates derived in this paper are new. Further, the duality argument in Sect. 5 seems to be the first of its kind that actually works without artificial assumptions on the regularity properties of dual quantities, cf. [31,33,43]. Compare also with the analysis and the counterexamples in [12] in this context, which demonstrate that the assumptions on the dual solution made in [31,33,40,43] are indeed unrealistic and cannot be expected to hold for the classical obstacle problem and which, in combination with our positive results in Sect. 5, show that it makes a huge difference for the behavior of the finite element error in the lower L p -norms whether the variational inequality under consideration involves inequality constraints on the boundary ∂Ω or in the interior of Ω.
To help the reader navigate this paper, we conclude this introduction with a brief overview of the structure and the content of the following sections: Section 2 is concerned with preliminaries. Here, we clarify the notation, state our precise assumptions, and collect several regularity results for the problem (1). In Sect. 3, we prove the Céa-type lemma and the supercloseness result that are at the heart of our error analysis. Section 4 addresses the consequences that the results of Sect. 3 have for the derivation of finite element error estimates. The main results of this section, Corollaries 4.1 and 4.2 and Theorem 4.3, contain various W 1, p (Ω)-, L p (Ω)-, and H 1/2 (∂Ω)-estimates that cover a large variety of different situations. Section 5 is devoted to the analysis of the L 4 -error in the case f ∈ L ∞ (Ω). Here, we extend the classical Aubin-Nitsche trick to (1) and prove that the continuous and the discrete solution satisfy u − u h L 4 (Ω) = O(h 2−ε ) for all ε > 0 when the involved contact sets are sufficiently well-behaved. In Sect. 6, we summarize our results and give some concluding remarks. Section 7 finally contains numerical experiments that confirm our theoretical findings.
With L k and H k , we denote the k-dimensional Lebesgue and Hausdorff measure, and with tr(·) the classical trace operator, cf. [3]. For functions v with a continuous representative, we typically drop the prefix tr and simply write v instead of tr(v). Further, we use the symbols cl(·) and ∂ to denote the topological closure and the boundary of a set, respectively, and the abbreviation B r (x) to denote the closed ball of radius r > 0 around an x ∈ R 2 . With O(·), o(·), we denote the classical Landau symbols, and with C a generic constant which may change within an estimate but is never dependent on crucial quantities as, e.g., the mesh width. If we want to emphasize that C depends on a quantity α, then we write C = C(α). Lastly, we define a + := max(0, a) and a − := min(0, a) for all a ∈ R.

The continuous setting
As already mentioned in the introduction, the purpose of this paper is to study finite element error estimates for the two-dimensional scalar Signorini problem Here and in what follows, Δ and ∂ n denote the (distributional) Laplacian and the normal derivative, respectively, and f ∈ L 2 (Ω) is a given right-hand side. To avoid obscuring the basic ideas of our analysis with technicalities and distinctions of cases and to reduce the notational overhead, throughout this paper, we always assume that Ω is the unit square, i.e., Ω := (0, 1) 2 . Our arguments can be extended straightforwardly to other convex polygonal domains with obvious modifications and, depending on the largest interior angle of Ω, possibly additional assumptions on the exponent p in the L p -and W 1, p -error estimates. The same is true for other variants of Signorini's problem as, e.g., the version studied in [2] whose partial differential operator does not contain a term of order zero and which involves separate Dirichlet-, Neumann-, and Signorini-boundary parts Γ D , Γ N , and Γ S . For such problems, however, some care has to be taken since the H 2regularity result in Proposition 2.1 and error estimates analogous to those in Lemma 3.2 are not directly available in the literature but first have to be established (under suitable assumptions on the angles between Γ D , Γ N , and Γ S ) by means of the techniques of [2,36,39]. We omit a detailed discussion of this topic to avoid overloading this paper. Lastly, we would like to point out that a Signorini problem of the type (2) with a sufficiently regular non-zero obstacle on the boundary can always be transformed into a problem with a vanishing obstacle by an elementary translation argument. This allows us to restrict our attention to the homogeneous situation in (2) without a great loss of generality.
To begin our analysis, we recall that the weak formulation of (2) is given by the elliptic variational inequality with the admissible set Note that Proposition 2.1 and the Sobolev embeddings imply that u possesses a representative which is continuous on the closure cl(Ω) of the domain Ω. In what follows, we always use this representative when talking, e.g., about level sets. As it turns out, solutions u to (S) enjoy additional regularity properties when the contact set {x ∈ ∂Ω | u(x) = 0} is sufficiently well-behaved. To explore this effect, we introduce: Definition 2.2 (Condition (A)) A solution u ∈ H 2 (Ω) of (S) is said to satisfy condition (A) if the relative boundary of the contact set {x ∈ ∂Ω | u(x) = 0} in ∂Ω has onedimensional Hausdorff measure zero and if the relative interior of the contact set in ∂Ω consists of at most finitely many connected components.
Under assumption (A), the solution u of the variational inequality (S) can be identified with the solution of a mixed Dirichlet-Neumann problem with segment-wise prescribed boundary conditions on a convex polygonal domain. This, together with the H 2 -regularity result in Proposition 2.1 and the analysis in [19], allows us to prove the following improved regularity result for the Signorini problem (S) which clarifies precisely which regularity assumptions on u are reasonable when it comes to the derivation of finite element error estimates: Theorem 2.3 (Improved regularity for Signorini's problem) Suppose that u ∈ H 2 (Ω) solves (S) and satisfies (A). Then, the following holds true: 1. If f is in L p (Ω) for some 2 < p < 4, then it holds u ∈ W 2, p (Ω).
Proof Since u satisfies condition (A), we may find relatively open disjoint straight line segments Γ i ⊂ ∂Ω, i = 1, . . . , N + M, N , M ∈ N 0 , and a set R ⊂ ∂Ω of one-dimensional Hausdorff measure zero such that From the variational inequality (S), the H 2 -regularity of the solution u, and Green's first identity, it follows further that Note that the parts of the contact set that are contained in the line segments Γ i , i = N + 1, . . . , N + M, are negligible here due to the properties of R. Let us suppose now that f is an element of L p (Ω) for some 2 < p < 4. Then, we may use [19,Theorem 4.4.3.7] to deduce that there exist real numbers c i and trigonometric functions φ i such that holds, where r i ≥ 0 and θ i ∈ [0, 2π) denote local polar coordinates centered at the vertices x i , i = 1, . . . , N + M, of the partition {cl(Γ i )} of the boundary ∂Ω, and where η i is a smooth cut-off function which is identical one in a neighborhood of x i for each i. We already know, however, that u ∈ H 2 (Ω), and it is easy to check that the factor r 1/2 i prevents a function of the form η i r . This implies that all c i in (5) have to be zero and proves the first claim, cf. also with the discussion in [41, Remark 2.1] and [32] in this regard. To obtain the second claim, we can proceed along exactly the same lines: If f is an element of L p (Ω) for some 4 < p < ∞, then we may use the same arguments as above and again [19,Theorem 4.4.3.7] to deduce that there exist real numbersc i and trigonometric functionsφ i with where η i , r i , and θ i are as in (5). The functions  Section 2] in this context, where it is supposed (but not explicitly stated) that each point in the relative boundary of the contact set admits an open neighborhood in which the solution u changes precisely once from contact to non-contact. (Note that the latter assumption implies in particular that the relative boundary of the contact set in ∂Ω is finite and is thus stronger than condition (A).) Theorem 2.3 is more precise in this regard in that it explicitly states which assumptions on the contact set are needed to rigorously prove improved regularity properties for (S) and further also quantifies which regularity can be expected for the "regular" part of the solution that remains when the singular contributions coming from the transition points on the boundary are subtracted from u. 3. For problems of the type (3) with right-hand side f ≡ 0 and an additional inhomogeneity on the boundary ∂Ω, it is possible to rigorously prove that condition (A) is satisfied. See the recent contribution [2] for details. 4. As already mentioned in the introduction, in the context of finite element error estimates for Signorini's problem, conditions similar to (A) have been needed for the derivation of optimal-order finite element error estimates in the H 1 -norm for a considerable amount of time. See, e.g., [5,8,9], for examples of contributions that require analogous assumptions on the contact set {x ∈ ∂Ω | u(x) = 0} of the continuous solution u. Only very recently in 2015, it was finally accomplished in [16] to get rid of these conditions and to establish H 1 -error estimates of optimal order for (S) that solely rely on standard Sobolev regularity properties. We would like to point out that, similar to the results in [16], the and H 1/2 (∂Ω)-error estimates that we derive in Sects. 3 and 4 for the problem (S) do not require condition (A), but only the regularity properties that are implied by it. Only the derivation of the L 4 -error estimates in Sect. 5 will make explicit use of condition (A) (and a comparable condition in the discrete setting). In what follows, in each theorem, lemma etc., we will clearly state whether the respective result requires condition (A) or suitable regularity properties of the solution u of the problem (S).

The discrete setting
As discrete counterparts of the variational inequality (S), we consider problems of the form Our standing assumptions on the quantities in (S h ) are as follows:

Assumption 2.5 (Standing assumptions for the FE-discretization)
See, e.g., [11,Definition 2] or [7,Definition 4.4.13] for the definition of the term "quasi-uniform family of triangulations". For brevity's sake, in what follows, we often ignore the upper bound on the mesh width and simply write "for all h > 0" instead of "for all 0 < h < 1/2". From standard results as, e.g., [25, Theorem II-2.1], we obtain: In the remainder of this paper, our aim will be to study the approximation properties of the solution u h of (S h ) for h 0. The main ingredients of our error analysis are:

A Céa-type lemma and a supercloseness result
To study the error u − u h , we introduce the following operator: Note that R h : H 1 (Ω) → V h is precisely the solution operator of the unconstrained problem associated with (S h ). In particular, R h is well-defined, linear, and continuous, and we may invoke classical results to obtain: with some constant C > 0 independent of h and v.
Proof For the domain Ω = (0, 1) 2 , the estimate (7) can be derived as follows: Consider an arbitrary but fixed v ∈ W 1,∞ (Ω) with associated Ritz projection R h (v) ∈ V h . Then, we may use reflections to extend v and R h (v) first to the rectangle (0, 1) × (−1, 2) and subsequently to the squareΩ : Here,R h is the Ritz projection operator on the mesh ofΩ that is obtained from the reflection procedure. From the interior norm estimate [38, Theorem 1.2], we may now deduce that there exists a Sinceṽ has the same W 1,∞ -norm as v, the above, the triangle inequality, and the properties ofṽ andR h (ṽ) imply that The above estimate, the regularity results of [19], and exactly the same arguments as in the proofs of [36, Equations (1.8), (1.9)] yield that, for every 2 It remains to prove the L p (∂Ω)-estimate. This, however, follows immediately from the last inequality and [19, Theorem 1.5.1.10] with parameter ε := h p . Note that the above argumentation only works for rectangles and squares. For more general domains Ω, (7) can be obtained by employing the techniques of [36,39]. To do so, however, one has to study in detail the behavior of certain Green's functions in the vicinity of the corners of the domain under consideration, cf. the comments in [39,Section 3]. Such a study is beyond the scope of this paper.
To prove the second assertion of the lemma, we suppose that a function v By proceeding completely analogously to the proof of [29,Lemma 5.7] (with the Besov estimate in [29,Lemma 4.1] replaced with the second estimate in [26, Lemma 2.1]), we obtain that, for every ε ∈ (0, 1/2), there exists a C > 0 with From the inequality of Poincaré-Wirtinger and the first part of the lemma, we may now deduce that holds with some C > 0. If we combine the above with (8), the trace theorem, the triangle inequality, and the L p (∂Ω)-estimate in (7), then the claim follows immediately.
We may now make the following observation (that has already been made in [11,Lemma 10] for the classical obstacle problem): This proves the claim.
The above result shows that it suffices to study the error that occurs in the solution u h of (S h ) when the original obstacle (i.e., the zero function) in (S h ) is replaced with the obstacle R h (u) − u to relate the approximate solution u h and the Ritz projection R h (u) of the exact solution u to each other. By pursuing this approach, we obtain the following Céa-type result: Lemma 3.4 (A Céa-type property) Let f ∈ L 2 (Ω) be arbitrary but fixed, and let u and u h denote the solutions of (S) and (S h ), respectively. Then, it holds Proof Consider an arbitrary but fixed w h ∈ V h that is contained in the set on the right-hand side of (9). (Note that this set is not empty since it contains R h (u).) Since On the other hand, we know that By addition, it now follows that This proves the claim.
Note that the above arguments work for all elliptic variational inequalities with unilateral constraints (not just for the Signorini problem). To obtain a tangible a priori estimate for the norm on ∂Ω and which has a small H 1 -norm. This is accomplished in the following lemma by means of a unilateral approximation technique that has also been used in [11,12,30,31,42]: is a function with a non-negative trace that can be decomposed into two parts v s and v r which satisfy the conditions in point 2 of Theorem 2.3 for some 4 < p < ∞. Then, for every ε ∈ (0, 1/2) and every Proof We first introduce some notation: Suppose that h > 0 is arbitrary but fixed. We denote the nodes of the triangulation T h which are contained in the boundary of the square Ω = (0, 1) 2 with x i , i = 0, . . . , N , starting with x 0 := (0, 0) and then proceeding counterclockwise. For convenience, we additionally set . . , N + 1, and τ i to be the mesh cell whose boundary contains σ i . With C we again denote a generic constant which may change within an estimate but never depends on h. We now proceed in three steps: Step 1 (h-Independent Morrey Inequality on the Mesh Cells) Consider the reference element T := conv{(0, 0), (1, 0), (0, 1)}. Then, we know from the classical Morrey inequality, see [1,Theorem 5.4], that for every 2 < p < ∞ there exists a constant From the inequality of Poincaré-Wirtinger, we obtain further that there exists another By combining the last two inequalities, we may deduce that holds for all v ∈ W 1, p (T ) with average value zero on T . Since the seminorms appearing here do not detect constant functions, the last estimate also holds for all v ∈ W 1, p (T ). Consider now an arbitrary but fixed Due to the quasi-uniformity of the underlying family of meshes, we can find a constant C independent of h and i with We may thus conclude that there exists a constant C > 0 independent of i and h with max Step 2 (Proof in the W 2, p -Case) Suppose now that a function v ∈ W 2, p (Ω), 2 < p < ∞, with a non-negative trace is given, let h > 0 be arbitrary but fixed, and consider the auxiliary problem where we use the C(cl(Ω))-representatives of v and R h (v). Since (11) is a finitedimensional minimization problem with a non-empty compact admissible set and a continuous objective functional, it admits at least one solutionw h ∈ V h . Consider now an arbitrary but fixed x i with 0 ≤ i ≤ N . Then, the fact thatw h solves (11) implies that we cannot reduce the function valuew h (x i ) (while leaving the other nodal values unchanged) without violating the constraint R h (v) − v ≤w h on ∂Ω. This implies that one of the following three cases has to hold true (as one may easily check by contradiction, cf. Figure 1 and the analysis in [11]): 3. There exists an a ∈ (x i , x i+1 ] satisfying (12).
In the second case, we may apply Taylor's formula in the direction of the line segment [ with some constant C independent of i and h. Here, we have used the properties of a, the regularity of v, the affine linearity ofw h and R h (v) on [x i−1 , x i ), and (10). Completely analogously, we obtain in the third case that We have now proved that and we may deduce by summation that Using the inverse estimate in [ with some constant C > 0 which may change from step to step but is always independent of h. To construct a function w h ∈ V h with the desired properties, it remains to extend tr(w h ) suitably to a function in V h . This can be accomplished, e.g., by employing the discrete harmonic extension operator E h : tr(V h ) → V h , which, according to [27,Lemma 3.2], satisfies for some constant C > 0 independent of h and v h .
Step 3 (Proof in the v s -v r -Case) For a function v with a non-negative trace that can be decomposed into two parts v s and v r which satisfy the conditions in point 2 of Theorem 2.3 for some 4 < p < ∞, we can proceed completely analogously to Step or one of the cases 2 and 3 at each node x i , i = 0, . . . , N . Let us again consider case 2, fix an ε ∈ (0, 1/2), write q := 2/(1 + 2ε) ∈ (1, 2), and assume w.l.o.g. that the line segment Then, we may use the same calculation as in (13), the regularity properties of v s and v r , and Morrey's inequality to obtain If we use exactly the same strategy in case 3, then it follows thatw h satisfies for all i = 0, . . . , N , where I h again denotes the Lagrange interpolation operator. By integration, we may now again deduce (using the estimate (a + b) q ≤ C(a q + b q )) and, by summation, with a constant C which may depend on ε but is independent of h. The claim now follows completely analogously to Step 2.
By combining Lemmas 3.4 and 3.5, we now arrive at the following main result of this section: Theorem 3.6 (Supercloseness) Suppose that u solves (S) for some f ∈ L 2 (Ω). Then, the following holds true for the Ritz projection R h (u) and the finite element solution u h of (S h ): 1. If u satisfies u ∈ W 2, p (Ω) for some 2 < p < 4, then there exists a constant C > 0 independent of h with 2. If u admits a decomposition u = u s + u r as in point 2 of Theorem 2.3 for some 4 < p < ∞, then, for every ε ∈ (0, 1/2), there exists a constant C > 0 independent of h with Note that Theorem 3.6 indeed shows that the Ritz projection R h (u) of the exact solution u of (S) is superclose to the finite element approximation u h characterized by (S h ) as, for the considered ranges of p, (14) and (15) yield estimates of the form with an exponent γ that is strictly greater than one. The H 1 -error between R h (u) and u h thus decays faster than that between u and u h which, in the considered situation, can be expected to be at most of order O(h). We would like to point out that the behavior, that we observe here, agrees very well with intuition. Since the additional constraint in (S h ) is only present on the boundary, it is only natural that the finite element solution u h is very close to the Ritz projection R h (u) which may be interpreted as the solution of an associated unconstrained problem.

Consequences for W 1,p -, L p -, and H 1/2 -error estimates
As the behavior of the error R h (u) − u is known by Lemma 3.2, Theorem 3.6 gives rise to estimates for the quantity u − u h in a straightforward manner. The results that are obtained along these lines are collected in the following two corollaries: Corollary 4.1 Suppose that u solves (S) for some f ∈ L 2 (Ω). Assume further that there exists a 2 < p < 4 with u ∈ W 2, p (Ω). Then, for every 1 < q < ∞, there exists a constant C > 0 independent of h with Proof From Theorem 3.6, the inverse estimates in [7, Theorem 4.5.11], Lemma 3.2, and the triangle inequality, it follows that This proves the first estimate. Similarly, we may compute (using standard error estimates for the Lagrange interpolant I h (u), see [7,Theorem 4.4.20], Sobolev embeddings, again [7, Theorem 4.5.11], and Lemma 3.2) that holds for all 1 < q < ∞. Further, we may use the discrete Sobolev inequality in [7, Lemma 4.9.2] and exactly the same arguments as above to obtain It remains to prove the H 1/2 (∂Ω)-error estimate. To this end, we define ψ h to be the unique element of V h , which is one at every boundary node and zero at every interior node, and supp(ψ h ) to be the support of ψ h . We may now use the classical trace theorem and Hölder's inequality to infer that Now we may proceed as before (using Lemma 3.
This proves the claim.

Corollary 4.2 Suppose that u solves (S) for some f ∈ L 2 (Ω). Assume further that u admits a decomposition u = u s +u r as in point 2 of Theorem 2.3 for some 4 < p < ∞.
Then, for all ε ∈ (0, 1/2), there exists a constant C > 0 independent of h with Proof The proof is completely along the lines of that of the last corollary and only requires some minor modifications. We include it for the sake of completeness: Note that the regularity properties of u s and u r imply that u ∈ W 2,q (Ω) holds for all q ∈ (2, 4), cf. [17,Theorem 6.5]. Consider now an arbitrary but fixed ε ∈ (0, 1/2). Then, we may invoke Theorem 3.6 and compute (using the same ideas as before) Analogously (by [7, Lemma 4.9.2]), (Note that the coefficient of ε in the exponent is completely unimportant here since we may always redefine ε.) Finally, we may compute (using the trace theorem, Lemma 3.2, and again (16)) This completes the proof.
The error estimates in Corollary 4.2 turn out to be of particular interest when the exponent p can be chosen arbitrarily large. Indeed, in this limit case, we obtain the following important result:

Theorem 4.3 (Optimal FE-estimates under regularity assumptions) Suppose that u solves (S) for some f ∈ L 2 (Ω). Assume further that u admits a decomposition u = u s + u r with the properties in point 2 of Theorem 2.3 for all 4 < p < ∞.
Then, for all ε ∈ (0, 1/2), there exists a constant C > 0 independent of h with Proof In the considered situation, we may apply Corollary 4.2 with an arbitrarily large p > 4 and Corollary 4.1 with a p which is arbitrarily close to four. The assertions of the theorem now follow immediately by invoking these results and by noting that, for all ε ∈ (0, 1/2), we have (due to the W 1,∞ -estimate and the W 1,q -estimate for all exponents 2 < q < 4) with some constant C independent of h. This completes the proof.
Several things are noteworthy regarding the last three results:  [16] in that they only rely on (realistic) assumptions on the Sobolev regularity properties of the exact solution u and do not require any additional conditions on the contact set {x ∈ ∂Ω | u(x) = 0}, cf. [8]. At least to the authors' best knowledge, for W 1, p -error estimates with p > 2, such results have not been available so far in the literature. The same seems to be the case for the L ∞ -and the W 1,∞ -error estimate in (17). Note that we have derived our L ∞ -error estimates without invoking the discrete maximum principle and without the associated assumptions on T h , cf. [11,13,34]. (17) has already been obtained in [41,Theorem 2.2] in dimensions two and three under the assumption that the solution u is in H 5/2−ε (Ω) for all ε ∈ (0, 1/2). Note that this regularity can only be expected to hold if the right-hand side f satisfies f ∈ H 1/2−ε (Ω) for all ε ∈ (0, 1/2). The regularity assumptions that we work with in our analysis differ from this in that they also allow for a W 2, p -part of the solution u and are thus realistic for general right-hand sides f ∈ L p (Ω), see Theorem 2.3. We further obtain the H 1/2 -estimate in (17) with arguments that seem to be more elementary than those in [41]. However, in contrast to the analysis in [41], our approach cannot be extended straightforwardly to the three-dimensional setting since unilateral approximation results analogous to those in Lemma 3.5 are only available in limited form in dimensions d ≥ 3, cf. the analysis in [11]. 5. In [41], an L 2 -error estimate of order 3/2 − ε is obtained as a corollary of the H 1/2 -error estimate on the boundary. We obtain this order of convergence even in the L ∞ -norm.

The H 1/2 -error estimate in
As the reader might have noticed, Theorem 4.3 only yields the suboptimal order of convergence 3/2 − ε for, e.g., the L 4 -error. We thus miss a factor h 1/2 in comparison with the approximation properties of the Lagrange interpolation operator. In what follows, we demonstrate that a better estimate can be obtained with a non-standard duality argument, and that the order two (minus epsilon), that is observed in numerical experiments, can also be recovered analytically provided the continuous solution u satisfies condition (A) and the contact sets of the finite element approximations u h behave sufficiently well in the limit h 0.

L 4 -error estimates of optimal order via an Aubin-Nitsche trick
To estimate the L 4 -error, we use an approach that has been proposed by Mosco in [31,Section 7] for the one-dimensional obstacle problem and consider two dual variational inequalities -one for each of the components (u − u h ) + = max(0, u − u h ) and (u − u h ) − = min(0, u − u h ).

A duality argument for the component (u − u h ) + under condition (A)
To formulate our first dual problem, we introduce the following notation: Note that, according to condition (A), at least for all sufficiently small h, the number of connected components of A • h is finite and equal to the number of connected components of A • . Given a u which solves (S) and satisfies (A) and a solution u h of (S h ), we now consider the following auxiliary problem: Here, Note that the solution z of (D) depends on h (since A • h and u h do). From standard results and the analysis in [19], we may deduce: Moreover, for all ε ∈ (0, 1/2) and all sufficiently small h > 0, z is in W 2,(4−ε)/3 (Ω) and satisfies with some constant C > 0 independent of h.
Proof The unique solvability of (D) for all h > 0 follows from [25, Theorem II-2.1]. Further, we may employ Stampacchia's lemma, see [3,Theorem 5.8.2], and use the test function v := z − ∈ L in (D) to deduce that This proves that we indeed have z ≤ 0 L 2 -a.e. in Ω and, as a consequence, that tr(z) = 0 holds H 1 -a.e. on A • h . Moreover, by choosing the test functions v = 0 and v = 2z in (D), and by exploiting the Sobolev embeddings, we obtain where C is the embedding constant of H 1 (Ω) → L (4−ε)/(1−ε) (Ω). This yields (18). It remains to prove the W 2,(4−ε)/3 -regularity of z and (19). To this end, we note that the non-positivity of z in Ω, the condition tr(z) ≥ 0 on A • h , and the properties of the set A • h imply that z is the unique weak solution of the problem Since A • h and its complement can be written as the union of at most finitely many straight line segments which meet at an angle π/2 or π , we may again invoke [19,Theorem 4.4.3.7] to deduce that z ∈ W 2,(4−ε)/3 (Ω) holds for all ε ∈ (0, 1/2). To obtain the estimate (19), let us assume that A • = ∅ and A • = ∂Ω (else the proof is trivial). In this case, condition (A) implies that A • consists of finitely many connected components, that the relative boundary of A • in ∂Ω consists of finitely many points b 1 , . . . , b N , N ∈ N, and that we may find a δ > 0 with dist(b i , b j ) > 4δ for all i = j and dist(b i , {(0, 0), (0, 1), (1, 0), (1, 1)}) > 4δ for all b i which are not themselves corner-points of the square Ω. Choose rotationally symmetric cut-off functions holds for all i = 1, . . . , N , where B r (b) denotes the closed ball of radius r > 0 around a b ∈ R 2 , and decompose z into the parts z 0 , z 1 , . . . , z N defined by Suppose further that h is so small that the set A • h \A • is contained in the union of the balls B δ (b i ), i = 1, . . . , N . (This is the case for all sufficiently small h due to the definition of A • h .) Then, δ and the functions ψ i are clearly independent of h, and we may compute that Here, we have used that A • ⊂ A • h , and that the rotational symmetry of the cut-off functions and the choice of δ imply ∂ n ψ i ≡ 0 on ∂Ω for all i. Note that the boundary conditions in the above problem are independent of h. We may thus invoke [19,Theorem 4.3.2.4] and (18) to deduce that, for every ε ∈ (0, 1/2), we have Here, C > 0 is a generic constant which depends on ε and ψ 0 but not on h and which may change from step to step.  2δ (a, 0), Here, we have again used the properties of z and the rotational symmetry of ψ i . Since z i and ψ i vanish outside of the ball B 2δ (a, 0), we may now deduce that the (trivial extension of) the functionz(x, y) := z i (x + a − τ h − 1/2, y) satisfies −Δz =ḡ L 2 -a.e. in Ω, If we expressz in terms of z i and use the same calculations as before, then we arrive at with a constant C > 0 which depends on ψ i and ε but is independent of h. Using the same arguments as above, we can transform each of the situations occurring at the points b i , i = 1, . . . , N , to one of finitely many reference configurations and use [19,Theorem 4.3.2.4] as well as (18) to prove that there exist constants C i independent of h with Note that, if we consider a point b i which is a corner of the square Ω, then the situation is even simpler than above since, in this case, the boundaries of A • and A • h are locally the same and equal to {b i } so that a translation argument as above is unnecessary. To arrive at (19), it now suffices to invoke the triangle inequality. This completes the proof.

Remark 5.3
Note that the solution z of the auxiliary problem (D) cannot be expected to possess W 2,q -regularity for some q ≥ 4/3 since it typically contains a singular part analogous to that in (5).
By choosing a suitable test function in (D) and by exploiting the estimates in Theorem 4.3, we may now deduce: Proposition 5.4 Suppose that u solves (S) for some f ∈ L ∞ (Ω) and that (A) is satisfied. Then, for all ε ∈ (0, 1/2), there exists a constant C > 0 independent of h such that, for all sufficiently small h > 0, we have Proof Let us denote the finitely many mesh nodes in the relative boundary of the set A • h with x i , i = 1, . . . , N , N ∈ N 0 , and the basis functions of the nodal basis of V h that belong to the nodes x i with ϕ i . Note that the number N is independent of h here for all sufficiently small h by our assumptions and the definition of where C is supposed to be a constant independent of h with diam(T ) ≤ Ch for all T ∈ T h . We claim that this v is admissible for (D). Indeed, on A • , we have u ≡ 0 and thus v ≥ 0. Further, we know that, for all small h and all x ∈ A • h \A • , we can find anx in the relative boundary of A • and a j ∈ {1, . . . , N } with x ∈ [x j ,x] and dist(x j ,x) < Ch, where [x j ,x] denotes the line segment between x j andx. This yields Choosing the function v + z ∈ L in (D) and using Lemma 5.2, the Sobolev embedding and Hölder's inequality now gives (with a generic C > 0 independent of h) Note that, since z vanishes on A • h , since the relative boundary of A • h consists of mesh nodes, and since z ≤ 0 in Ω, the Lagrange interpolant I h (z) ∈ V h satisfies I h (z) = 0 on A • ⊂ A • h and I h (z) ≤ 0 in Ω. This implies in combination with (S h ) and the reformulation (4) of (S) that i.e., we have Using the last inequality, standard results for the Lagrange interpolation operator, and again Lemma 5.2, we can continue the estimate in (20) as follows The above yields, in combination with Theorems 4.3 and 2.3, that there exists a constant C independent of h with where the Landau symbol refers to the limit ε 0. Using the above in (21) and performing the same calculation as before yields This proves the claim (after redefining ε).

A duality argument for the component (u − u h ) −
To obtain an L 4 -error estimate for the component (u − u h ) − , we can proceed along roughly the same lines as in the last subsection provided the contact sets of the finite element solutions u h behave sufficiently well for h 0.
To be more precise, we need the following assumption: Analogously to the last section, we may now consider the following auxiliary problem:z By invoking again the results of [19], we obtain: Moreover, for all ε ∈ (0, 1/2) and all sufficiently small h > 0,z is in W 2,(4−ε)/3 (Ω) and satisfies with some constant C > 0 independent of h.
Proof The unique solvability of (D) for all h > 0 follows from [25, Theorem II-2.1], and the non-positivity ofz and the property tr(z) = 0 onÃ h are obtained completely analogously to the proof of Lemma 5.2. The same is the case for the estimate (22). It remains to prove the W 2,(4−ε)/3 -regularity ofz and (23) for all sufficiently small h > 0. The former follows immediately from (A h ), [19,Theorem 4.4.3.7] and the same arguments as in Lemma 5.2. To obtain the latter, we assume that h is so small that the conditions in (A h ) hold with some d i ∈ ∂Ω, δ i > 0, i = 1, . . . , N , N ∈ N (for N = 0 the claim is trivial) and choose rotationally symmetric cut-off functions for all i, and such that the sets supp(ψ i ) ∩ ∂Ω have non-zero distance from each other and the corners of the square Ω. In this situation, the properties of the functions ψ i imply that we may find another cut-off function φ ∈ C ∞ c (R 2 ) with 0 ≤ φ ≤ 1 in R 2 and φ ≡ 1 in a neighborhood of the boundary ∂Ω such that the supports of the functionsψ i := ψ i φ are disjoint. Using theseψ i , we decomposez into the partsz 0 ,z 1 , . . . ,z N defined bỹ Since (A h ) implies that the relative boundary ofÃ h is contained in the union of the balls B δ i (d i ), i = 1, . . . , N , that each set B δ i (d i ) ∩ ∂Ω contains precisely one point of the relative boundary ofÃ h , and that the connected components ofÃ h each have a non-empty relative interior, we may argue as in the proof of Lemma 5.2 to deduce thatz 0 satisfies with some closed set B ⊂ ∂Ω whose connected components each have a nonempty relative interior and whose relative boundary consists precisely of the points d 1 , . . . , d N . Note that it is (at least in theory) possible that B varies with h since it is not uniquely determined by the above conditions. However, it is easy to check that only two sets B and combinations of boundary conditions are possible here. We may thus again invoke [19,Theorem 4.3.2.4] and use (22) to deduce that, for every ε ∈ (0, 1/2), we have with a generic constant C > 0 which is independent of h. It remains to estimate the W 2,(4−ε)/3 -norm of the functionsz i , i = 1, . . . , N . So let us consider an arbitrary but fixed point d i . Since d i is not a corner of Ω by (A h ), we may assume w.l.o.g. that d i = (a, 0) holds for some a ∈ (0, 1). Further, it follows from a straightforward calculation and the properties ofψ i that in Ω, Since the support supp(ψ i ) contains precisely one point (ã, 0) ∈ ∂Ω of the relative boundary ofÃ h by our assumption (A h ), we may complement (24) with one of the boundary conditions Using exactly the same arguments as in the proof of Lemma 5.2, we may now transform the situation at d i into one of finitely many reference configurations and invoke [19,Theorem 4.3.2.4] as well as (22) to deduce that there exists a constant C i independent of h with Proceeding exactly along the same lines at the other points d i and using the triangle inequality, we arrive at (23). This completes the proof.
By choosing the test function v =z + u − u h ∈L in (D), we now obtain: Proof Note that the definition ofÃ h implies u − u h = u ≥ 0 onÃ h . The function u − u h is thus an element ofL and we may choose the function v =z + u − u h in (D) to obtain Since the setÃ h consists of cells of the boundary mesh, sincez vanishes inÃ h , and sincez ≤ 0 holds L 2 -a.e. in Ω, we know that the Lagrange interpolant I h (z) vanishes inÃ h and that I h (z) is non-positive everywhere. This implies that, for all small enough s > 0, we have u h ± s I h (z) ∈ K h and, as a consequence, that Using the above in (26) yields for all ε ∈ (0, 1/2). A calculation completely analogous to that at the end of the proof of Proposition 5.4 now yields (25) as claimed.

Summary of results and remarks on the error analysis
In summary, we have now proved the following for problems (S) whose right-hand sides are in L ∞ (Ω) and whose solutions satisfy condition (A): Theorem 6.1 (Optimal FE-estimates under assumption (A)) Suppose that u solves (S) for some f ∈ L ∞ (Ω) and that condition (A) is satisfied. Then, for all ε ∈ (0, 1/2), there exists a constant C > 0 independent of h with If, additionally, the approximate solutions u h satisfy (A h ), then we also have  (23). Showing that the dual solution possesses W 2,(4−ε)/3 -regularity is relatively simple. 5. Recall that, for a classical obstacle problem with an essentially bounded righthand side, it can be shown that the exact solution enjoys W 2, p -regularity for all 2 ≤ p < ∞, cf. [11,25,28]. This implies in particular that it is possible to prove L ∞ -error estimates of order O(h 2−ε ) for arbitrarily small ε > 0. For the Signorini problem, this is different since the exact solution u cannot be expected to be in W 2,4 (Ω) even for smooth right-hand sides f . The L 4 -error estimates in Theorem 6.1 thus yield an order of convergence that cannot be recovered with pointwise a priori error estimates. 6. Note that Theorem 6.1 shows that a counterexample, which demonstrates that the L 4 -error u − u h L 4 (Ω) is in general not of order O(h 2−ε ), has to be very exotic (if it exists) since the contact sets of u and u h have to exhibit a very degenerate behavior for the conditions (A) and (A h ) to be violated.

Numerical experiments
We conclude this paper with numerical experiments that confirm our theoretical findings. To construct a model problem that allows us to validate our results and that possesses a known analytic solution, we proceed along the lines of [41, Section 7] and consider the functionũ : R × (0, ∞) → R, x → −r 3/2 sin( 3 2 θ). Here, r and θ denote polar coordinates centered at (0.5, 0), i.e., Note that the functionũ is exactly of the same type as the singular terms on the righthand side of (6), cf. also the analysis in [19,20]. Moreover, it is easy to check that u is an element of (In the experiments below, this ψ was an appropriately defined ninth-order spline.) Then, the properties ofũ and ψ yield that the map u : Ω → R, x → 10 ψ(r )ũ(r , θ), where Ω again denotes the unit square (0, 1) 2 and where we have added the factor ten for scaling reasons. Note that the conditions in (29) imply in particular that the function u solves (S) with right-hand side f := −Δu + u ∈ C(cl(Ω)). What we have constructed in (28) is thus indeed an analytic solution of Signorini's problem that can be used as a reference in our numerical experiments, cf. Fig. 2.
The results that we have obtained for the right-hand side f associated with the solution u in (28) by means of the finite element scheme described in Sect. 2.3 can be seen in Tables 1 and 2. Here, we have used Friedrichs-Keller triangulations to discretize the continuous problem (S) and a 16-point Gauss-Legendre-type quadrature rule for triangles to evaluate the various W s, p -errors and the integrals arising on the right-hand side of (S h ). The finite-dimensional elliptic variational inequalities obtained from the discretization scheme have been solved by means of the MatlabR2019b-routine quadprog to a high precision. Note that the choice of mesh widths in our numerical The H 1/2 -error was estimated by means of the interpolation property of the space H 1/2 (∂Ω) The row "regr." contains the EOCs that are obtained from the data in Table 2 by linear regression (after a logarithmic scaling) and the row "theo." the orders of convergence expected from our analysis. The expected rates of convergence for the L 16 -and the W 1,16 -error have been computed with inverse estimates experiments ensures that the point (0.5, 0), where the analytic solution u detaches from the boundary ∂Ω and where ∇ 2 u possesses a singularity, never coincides with a node of the underlying mesh. This constitutes the worst case scenario as the critical part of the contact set of u is never resolved properly. Further, it should be noted that the condition (A) is trivially satisfied in (28) by the properties of the analytic solution u. Our numerical experiments indicate that the same is true for the condition (A h ) as it can be observed that the contact sets of the finite element solutions u h approximate their continuous counterpart, cf. the comments after Definition 5.5. As the results in Tables 1 and 2 illustrate, the behavior observed in our numerical experiments agrees very well with the predictions in Theorem 6.1. (Note that this result is indeed applicable here since f ∈ C(cl(Ω)).) In particular, the experimental orders of convergence (EOCs), i.e., the quantities fit very well to the a priori error estimates in (27). Table 2 further shows that the rates of convergence in the L p -and the W 1, p -norms break down in the situation of (28) when p is greater than the critical exponent four. This demonstrates that, for instance, the order one (minus epsilon) is in general unobtainable when we consider the W 1, p -error for some p > 4. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.