Existence, iteration procedures and directional differentiability for parabolic QVIs

We study parabolic quasi-variational inequalities (QVIs) of obstacle type. Under appropriate assumptions on the obstacle mapping, we prove the existence of solutions of such QVIs by two methods: one by time discretisation through elliptic QVIs and the second by iteration through parabolic variational inequalities (VIs). Using these results, we show the directional differentiability (in a certain sense) of the solution map which takes the source term of a parabolic QVI into the set of solutions, and we relate this result to the contingent derivative of the aforementioned map. We finish with an example where the obstacle mapping is given by the inverse of a parabolic differential operator.


Introduction
Quasi-variational inequalities (QVIs) are versatile models that are used to describe many different phenomena from fields as varied as physics, biology, finance and economics. QVIs were first formulated by Bensoussan and Lions [8,31] in the context of stochastic impulse control problems and since then have appeared in many models where nonsmoothness and nonconvexity are present, including superconductivity [28,4,36,7,39], the formation and growth of sandpiles [38,5,37,35,7] and networks of lakes and rivers [37,35,6], generalised Nash equilibrium problems [23,17,33], and more recently in thermoforming [2].
In general, QVIs are more complicated than variational inequalities (VIs) because solutions are sought in a constraint set which depends on the solution itself. This is an additional source of nonlinearity and nonsmoothness and creates considerable issues in the analysis and development of solution algorithms to QVI and the associated optimal control problems.
We focus in this work on parabolic QVIs with constraint sets of obstacle type and we address the issues of existence of solutions and directional differentiability for the solution map taking the source term of the QVI into the set of solutions. This extends to the parabolic setting our previous work [2] where we provided a differentiability result for solution mappings associated to elliptic QVIs. Literature for evolutionary QVIs is scarce in the infinite-dimensional setting: among the few works available, we refer to [24] for a study of parabolic QVIs with gradient constraints, [25] for existence and numerical results in the hyperbolic setting, [18] for QVIs arising in hydraulics, [27] for state-dependent sweeping processes, and evolutionary superconductivity models in [39], as well as the work [19]. Differentiability analysis for parabolic (non-quasi) VIs was studied in [26] and [15].
Let us now enter into the specifics of our setting. Let V ⊂ H ⊂ V * be a Gelfand triple of separable Hilbert spaces with V c ֒− → H a compact embedding. Furthermore, we assume that there exists an ordering to elements of H via a closed convex cone H + satisfying H + = {h ∈ H : (h, g) ≥ 0 ∀g ∈ H + }; the ordering then is ψ 1 ≤ ψ 2 if and only if ψ 2 − ψ 1 ∈ H + . This also induces an ordering for V and V * as well as for the associated Bochner spaces L 2 (0, T ; H), L 2 (0, T ; V ) and so on. We define V + := {v ∈ V : v ≥ 0}. We write v + for the orthogonal projection of v ∈ H onto the space H + and we have the decomposition v = v + − v − . Suppose that v ∈ V implies that v + ∈ V and that there exists a constant C > 0 such that for all v ∈ V , v + V ≤ C v V . An example of such a space V is V = W 1,p (Ω) for 1 ≤ p ≤ ∞; we refer to [1] for a definition of the Sobolev space W 1,p (Ω) over a domain Ω.
Further assumptions will be introduced later as and when required. We consider parabolic QVIs of the following form.
QVI Problem: Given f ∈ L 2 (0, T ; H) and z 0 ∈ V , find z ∈ L 2 (0, T ; V ) with z ′ ∈ L 2 (0, T ; V * ) such that We write the solution mapping taking the source term into the (weak or strong) solution as P z0 so that (1) reads z ∈ P z0 (f ) (sometimes we omit the subscript). In this paper, we contribute three main results associated to this QVI: • Existence of solutions to (1) via time-discretisation (Theorem 2.9): we show that solutions to (1) can be formulated as the limit of a sequence constructed from considering time-discretised elliptic QVI problems. This result makes use of the theory of sub-and supersolutions and the Tartar-Birkhoff fixed point method.
• Approximation of solution to (1) by solutions of parabolic VIs (Theorem 3.8): we define an iterative sequence of solutions of parabolic VIs and show that the sequence converges in a monotone fashion to a solution of the parabolic QVI; this is another QVI existence result. The existence for the aforementioned parabolic VIs comes from either Theorem 2.9 or from a certain result of Brezis (which will be given in the relevant section), giving rise to two different sets of assumptions under which the theorem holds.
• Directional differentiability for the source-to-solution mapping P (Theorem 5.15): we prove that the map P is directionally differentiable in a certain sense using Theorem 3.8 and some technical lemmas related to the expansion formulae for parabolic VI solution mappings.
Thus we give two existence results and a differential sensitivity result. It should be noted that the differentiability result essentially gives a characterisation of the contingent derivative (a concept frequently used in set-valued analysis) of P (between appropriate spaces) in terms of a parabolic QVI; see Proposition 5.16 for details.

Notations and layout of paper
We shall usually write the duality pairing between V and V * as ·, · rather than ·, ·, V * ,V for ease of reading. We use the notations ֒→ and Note the relationships W s (0, T ) ֒→ W (0, T ) and W r (0, T ) ⊂ W (0, T ). If z ∈ W (0, T ) satisfies (1), we say that it is a weak solution or simply a solution. A weaker notion of solution is given by transferring the time regularity of the solution onto the test function and requiring only z ∈ L 2 (0, T ; V )∩L ∞ (0, T ; H) to satisfy and we call z a very weak solution. Note that the initial data also has been transferred onto the test function (indeed, the weak solution is not sufficiently regular to have a prescribed initial data).
The paper is structured as follows. In §2 we consider the existence of solutions to (1) via the method of time discretisation: we characterise a solution of the parabolic QVI as the limit of solutions of elliptic QVIs. In §3, we approximate solutions of the QVI by a sequence of solutions of parabolic VIs that are defined iteratively. In §4 we consider parabolic VIs and directional differentiability with respect to perturbations in the obstacle and we make use of this in §5 to prove that the source-to-solution map P is directionally differentiable in a particular sense. We highlight some possible alternative approaches in §6 and finish with an example in §7.

Existence for parabolic QVIs through time discretisation
We prove existence to (1) by the method of time discretisation via elliptic QVIs of obstacle type. This is in contrast to [25] where the discretisation for evolutionary QVI was done in such a way as to yield a sequence of elliptic VIs, and as far as we are aware, our approach is novel in these type of problems.
We make the following basic assumption.
Let N ∈ N, h N := T /N and for n = 0, 1, ..., N , t N n := nh N . This divides the interval [0, T ] into N subintervals of length h N ; we will usually write h for h N . We approximate the source term by and we consider the following elliptic QVI problem.
Discretised problem: Given z N 0 := z 0 , find z N n ∈ V such that This problem is sensible since by (3), for fixed t, the mapping v → Φ(v)(t) is well defined from V into V , since we may consider V ⊂ L 2 (0, T ; V ) (elements of V can be thought of constant-in-time elements of L 2 (0, T ; V )), and Φ(v) can be evaluated pointwise in time. We consider first the existence of solutions to (4).

Existence and uniform estimates for the elliptic approximations
The inequality (4) can be rewritten as We write the solution of (4) or (5) as Q t N n−1 (hf N n + z N n−1 ) ∋ z N n . Related to (5) is the following VI: denote by E t,h (g, ψ) = v the solution of this problem. For fixed t and h and sufficiently smooth data, this mapping is well defined since the obstacle Φ(ψ)(t) ∈ V . To ease notation, we write E t N n (g, ψ) instead of E t N n ,h N (g, ψ) (where h N is the step size) because specifying h N is redundant.
Let us make an observation which follows from the theory of Birkhoff-Tartar [41,12]. Suppose there exists a subsolution z sub and a supersolution z sup to (4), i.e., z sub ≤ E t N n−1 (hf N n + z N n−1 , z sub ) and z sup ≥ E t N n−1 (hf N n + z N n−1 , z sup ). Then the QVI problem (4) has a solution z N n = E t N n−1 (hf N n + z N n−1 , z N n ) ∈ [z sub , z sup ]. The next lemma applies this idea. First, let us set A h w := w + hAw (this A h is the elliptic operator appearing in (5)) and definez n,N as the solution of the following elliptic PDE: Then the approximate problem (4) has a non-negative solution z N n ∈ [z N n−1 ,z n,N ]. Thus the sequence {z N n } n∈N is increasing.
Proof. We first show that z N 0 = z 0 is a subsolution for the QVI for z N 1 . Indeed, let w : again by assumption (D2). This shows that . Then we apply the theorem of Birkhoff-Tartar to obtain with the final inequality because the obstacle associated to t N n is greater than or equal to the obstacle associated to t N n−1 by assumption (6). We also havez are sub-and super-solutions respectively for Q t N n (f N n+1 + z N n ) (and the supersolution is greater than the subsolution). Therefore, by Birkhoff-Tartar there exists a z N n+1 It appears that we may select any solution z N n as given by the Tartar-Birkhoff theorem in the previous lemma, regardless of how we choose z N n−1 . For example, we may choose z N n−1 to be the minimal solution on its corresponding interval (with endpoints given by the sub-and supersolution) and z N n to be the maximal solution on its corresponding interval, with no effect on the resulting analysis (though of course, different choices may lead to different solutions of the original parabolic QVI in question in the end).
We now obtain some bounds on the sequence {z N n }. For this, we use the fact that if f ∈ L 2 (0, T ; H), then to get Neglecting the second term and summing this up from n = 1 to n = m and using (7), we obtain

Interpolants
We define the piecewise constant interpolants where χ A represents the characteristic function on the set A. In order to ease the presentation, we often use the notation Proof. Since the T N n are disjoint, we see that (8), and (9). A similar argument leads to the bounds on z N − .
To be able to handle the time derivative, it is useful to construct the interpolant known as Rothe's function, which also has the time derivative Corollary 2.5. Under the assumptions of Lemma 2.3 Proof. We see that with the last two inequalities by (9). Regarding the time derivative, using (10), we find

Passing to the limit in the interpolants
By the previous subsection, we have the existence of z,ẑ such that, for a relabelled subsequence, the following convergences hold: Lemma 2.6. We haveẑ ≡ z. Furthermore, z N − ⇀ z in L 2 (0, T ; V ) and weakly-star in L ∞ (0, T ; H).

Proof. Observe that
with the final line by (10). This shows that z N − −ẑ N → 0 in L 2 (0, T ; H), allowing us to identifyẑ = z as desired.
The convergence results above obviously imply thatẑ N → z in C 0 ([0, T ]; H) (due to Aubin-Lions), so that z 0 = z N (0) → z(0), i.e., z has the right initial data. Let us now see that z is feasible.
Lemma 2.7 (Feasibility of the limit). Let the following conditions hold: Then Proof. Since by (6), Using the assumption (12) applied to the right-hand side above, we have Passing to the limit here using (13) gives the result.
Regarding the assumptions of the previous lemma, (13) is a mild continuity requirement on Φ whereas for (12), and now supposing t ∈ T N j−1 for some j, this becomes and since j is arbitrary, this holds for all t. Hence (12) holds with equality. In order to pass to the limit in the inequality satisfied by the interpolant z N , we have to be able to approximate test functions in the limiting constraint set. This is possible as the next lemma shows.
Then for every v ∈ L 2 (0, . by assumption (14) as long as N ≥ N 0 . This shows that v N → v.
Let us consider two cases under which the assumption (14) of the previous lemma holds. (14) translates to a complete continuity assumption. Indeed, the sum term in (14) is and so (14) also holds.
Proof. Let v ∈ L 2 (0, T ; V ) satisfy v(t) ≤ Φ(z)(t) and let us take v N as defined in the proof of Lemma 2.8. Then by (4), Writing the duality product involving the time derivative as an inner product, we have, using the convergences in (11) and the weak lower semicontinuity of the bilinear form generated by A, so that z ∈ P(f ) ∩ W s (0, T ). Since {z N n } are non-negative, it follows that z is too. By (11), it follows that z Nj (t) → z(t) in H for almost every t ∈ [0, T ]. Let now s ≤ r and suppose that s ∈ T  i } i∈N is increasing. Passing to the limit, we find for almost every r and s with s ≤ r that z(s) ≤ z(r), i.e., t → z(t) is increasing.

Parabolic VI iterations of parabolic QVIs
In this section, we will show the existence of sequences of solutions to VIs that converge to solutions of QVIs. We begin with collecting some facts regarding parabolic VIs.
Consider the parabolic VI We write the solution as z := σ z0 (f, ψ) when it exists. Given f ∈ L 2 (0, T ; H) and ψ ∈ V independent of time, the solution z ∈ W s (0, T ) exists uniquely, see e.g. [3,11,9]. The problem (15) can be transformed to a parabolic VI with zero initial data with right-hand side f −Az 0 and obstacle ψ − z 0 with the substitution u(t) = z(t) − z 0 : We often write simply σ rather than σ z0 when we do not need to emphasise the initial data. The next lemma shows that σ is increasing in its arguments.
Proof. The z i satisfy the inequalities Let us take here v 1 = z 1 − (z 1 − z 2 ) + , which is clearly a feasible test function, and take v

This gives us
Adding yields whence using T-monotonicity and coercivity, we obtain z 1 ≤ z 2 .
We start by giving some existence results for (15). The first one is a result of applying Theorem 2.9 of §2 using the obstacle mapping Φ(v)(t) ≡ ψ(t) (it can be seen that all of the assumptions of the theorem are satisfied, refer also to the remarks below the proof of Lemma 2.8). The next proposition applies a result due to Brezis-Stampacchia (see [29, §2.9.6.1, p. 286]) to obtain existence of a very weak solution and then a further argument is required to obtain additional regularity. Proposition 3.3 (Existence II). Let ψ ∈ W r (0, T ) be such that t → ψ(t) is increasing with z 0 ≤ ψ(0). Then (15) has a unique solution z = σ(f, ψ) ∈ W r (0, T ).
which, by standard parabolic theory, has a strong solution v ǫ ∈ W s (0, T ) thanks to the regularity on ψ. A rearrangement and adding and subtracting the same term leads to Testing the equation above with (v ǫ − ψ) + and using the non-negativity of the right-hand side, Then [40,§III,Proposition 7.2] implies that the solution is actually strong, i.e., (15) holds and it belongs to W (0, T ) with the additional regularity Lz ∈ L 2 (0, T ; H), i.e., z ∈ W r (0, T ).
Some related regularity results given sufficient smoothness for f can be found in e.g. [10, Theorem 2.1].

Proposition 3.4 (Existence for
and either Then (16) has a solution z = S(f, ψ) ∈ W (0, T ) with the regularity that, in the first case, z ∈ W s (0, T ) is non-negative and t → z(t) is increasing, whereas in the second case, z ∈ W r (0, T ).
Proof. The first set of assumptions imply that the hypotheses of Proposition 3.2 hold for the obstacle Φ(ψ), whilst under the second set of assumptions, we apply Proposition 3.3.
The next lemma states that the solution mapping S(f, ψ) = z is increasing with respect to the arguments. This follows simply by using Lemma 3.1 and the fact that Φ is increasing.

Iteration scheme to approximate a solution of the parabolic QVI
We say that a function z sub ∈ L 2 (0, T ; V ) is a subsolution for (1) if z sub ≤ S(f, z sub ), and a supersolution is defined with the opposite inequality.
Remark 3.6. Let Φ(0) ≥ 0. If f ≥ 0, the function 0 is a subsolution, and the functionz defined bȳ is a supersolution to (1). Both claims follow by the comparison principle: the first claim is clear and the second follows upon realising thatz = S(f, ∞) ≥ S(f,z). We need the sign condition on f for z sub = 0 ≤ z sup =z.
Proof. This is obvious The previous lemma tells us in particular that any element of P(f ) is a subsolution for P(f + sd). The next theorem, which shows that a solution of the QVI can be approximated by solutions of VIs defined iteratively with respect to the obstacle, is based on an iteration idea of Bensoussan and Lions in [11, Chapter 5.1] (there, the authors consider Φ to be of impulse control type). The theorem and its sister result Theorem 3.10 show in particular that the approximating sequences converge to extremal (the smallest or largest) solutions of the QVI in certain intervals.

Theorem 3.8 (Increasing approximation of the minimal solution of QVI by solutions of VIs
Let either , hold. Then the sequence {z n } n∈N denoted by is well defined, monotonically increasing and satisfies

Furthermore,
• in the first case, z n ∈ W s (0, T ) with z n ≥ 0 , ∂ t z n ⇀ ∂ t z in L 2 (0, T ; H) and z ∈ W s (0, T ) is a strong solution (i.e. it satisfies (1) with additional regularity on the time derivative), and both z n and z are increasing in time • or, in the second case, z n ∈ W r (0, T ).
Before we proceed, let us observe that Proof. The proof is split into five steps. Lemma 3.5). This shows that z n is a monotonically increasing sequence.
2. EXISTENCE OF {z n }. For the actual existence, we apply Proposition 3.4 as we see now.
Let us also see why , which along with (24) implies that we can take the trace of the previous inequality where the first inequality is with the aid of (18). Making use of (26), the increasing property and (25) (which tells us that t → Φ(z 1 )(t) is increasing, since t → z 1 (t) is), Proposition 3.4 is again applicable and we use it to obtain z 2 .
By bootstrapping this argument we get that z n ∈ W s (0, T ) is well defined. Furthermore we have z n ≥ 0 by the sign condition on the data.
Second case. In case of (27) and (28), (18), (17) and (21) and (18) hold for the obstacle z sub and we get Suppose that the first part of (29) holds. Then since z 1 (0) = z 0 , we find that z 0 ≤ Φ(z 1 )(0) and then again (D2) is satisfied for the obstacle Φ(z 1 ), giving existence of z 2 = S(f, z 1 ). Repeating this, we get z n ∈ W r (0, T ). If instead the second part of (29) holds, we get by monotonicity that z 0 ≤ z 1 and using the increasing property of Φ, z 0 ≤ Φ(z 0 )(0) ≤ Φ(z 1 )(0) (the first inequality by (22)) and again we can apply the existence and proceed in this manner for general n.
3. UNIFORM BOUNDS ON {z n }. By (O1a) and the fact that z n ≥ 0, or by (O1b), we find that 0 is a valid test function in (30) and testing with it yields which immediately leads to bounds in L ∞ (0, T ; H) and L 2 (0, T ; V ) giving the weak convergences stated in the theorem to some z, for the full sequence (and not a subsequence) thanks to the monotonicity property.

Uniform bounds on ∂ t z n under first set of assumptions.
In this case, we can obtain a bound on the time derivative. Indeed, due to the work on the time discretisations in §2, from (11) and Lemma 2.6 we know that for each j, the interpolant One should bear in mind that the {(z j ) N } N are interpolants constructed from solutions of elliptic VIs and not QVIs since the {z j } j∈N are solutions of VIs. One observes the bound where the constant C ultimately arises from (10). Evidently, it only depends on the initial data and the source term. Hence, in this case, 4.1. Under first set of assumptions. In this case, using the strong convergence v n → v * assured by (O2a), we can use (31) and pass to the limit after writing the duality pairing for the time derivative as an inner product to get the inequality

Under second set of assumptions.
In this case, we take the limiting test function v * ∈ W (0, T ) with v * (0) = z 0 and rewrite (30), using the monotonicity of the time derivative and assumption (O3b) (which guarantees that Φ(z) ∈ W (0, T ), and hence v n ∈ W (0, T )) as By (O3b), we find that v n → v * in W (0, T ), and hence we can pass to the limit in the above to obtain (2).

MINIMALITY OF THE SOLUTION.
Suppose that z * ∈ P(f ) is the minimal solution on the interval [z sub , ∞), which in particular implies z * ≤ z. We see by the comparison principle and since z 0 = z sub is a subsolution that By this, we find Passing to the limit shows that z ≤ z * so that z = z * .
Remark 3.9. Note that the compactness assumptions (O2a) and (O3b) are only required for identifying the limit point z and showing that it is feasible.
Theorem 3.10 (Decreasing approximation of the maximal solution of QVI by solutions of VIs). Let z 0 := z sup be a supersolution of P(f ) and assume that Under the assumptions of the previous theorem (except (29)) except with z sub replaced with z sup and (27) replaced with the sequence {z n } is monotonically decreasing and converges to a solution z ∈ P(f ) with the same regularity and convergence results as stated in Theorem 3.8. Furthermore, z is the maximal solution of (1) or (2) in the interval Proof. It follows that z ∈ P(f ) ∩ (−∞, z sup ] by the same argumentation as in the proof of the previous theorem. Let us prove the claim of the maximality of the solution. Suppose that there exists a maximal solution z * ∈ (−∞, z sup ] so that z * ≥ z where z = lim n z n with z 0 = z sup . We have z 0 = z sup ≥ S(f, z sup ) ≥ S(f, z * ) = z * , and thus whence passing to the limit, z ≥ z * , and thus z = z * is the maximal solution.

Transformation of VIs with obstacle to zero obstacle VIs
It will become useful to relate solutions of the parabolic VI (16) with non-trivial obstacle to solutions of parabolic VIs with zero (lower) obstacle. We achieve this as follows. Take w 0 ≥ 0 and defineS w0 : w the solution to the parabolic VI with lower obstacle Omitting when convenient the subscript, we obtain the following estimate for w i =S(g i ): which then leads to and (due to Young's inequality applied to the right-hand side) Let us note that (32) in particular implies The relationship between solutions of VIs with non-trivial obstacles and VIs with zero obstacle is given in the next result. (19), and (20) hold, then in fact S z0 (g, ψ) ∈ W s (0, T ) and hence the spatial regularity AS z0 (g, ψ) ∈ L 2 (0, T ; H).

Expansion formula for variations in the obstacle and source term
The aim in this section is obtain differential expansion formulae for the solution mapping of parabolic VIs with respect to perturbations on the source term and the obstacle. This will form the backbone of our QVI differentiability result in the next section.

Definitions and cones from variational analysis
To state directional differentiability results for VIs, we need some concepts and notation which we shall collect in this subsection. Let us define the lower obstacle sets and for y ∈ K 0 , the radial cone at y We shall consider K 0 as a subset of the Banach space L 2 (0, T ; V ). The tangent cone is defined as the closure of the radial cone: . We now show that the tangent cone is contained in a set which has a convenient description (see also the discussion after Remark 5.6 in [15]).

whereȳ(t) is a quasi-continuous representative of y(t).
Proof. If w ∈ T rad K0 (y), then from (36), w ∈ L 2 (0, T ; V ) and there exists ρ * ≥ 0 such that y(t) + ρw(t) ∈ K 0 for almost every t and for all ρ ∈ [0, ρ * ), meaning that w(t) ∈ T rad K0 (y(t)) ⊂ T tan K0 (y(t)) for almost every t. This shows that and if we take the closure in L 2 (0, T ; V ) on both sides, Suppose that {w n } ⊂ L 2 (0, T ; V ) is a sequence that belongs to the set on the right-hand side above with w n → w in L 2 (0, T ; V ). Thus, for a subsequence, w nj (t) → w(t) in V and w nj (t) ∈ T tan K0 (y(t)) for almost every t. Since the tangent cones are closed sets, the limit point T ]} and the closure can be omitted on the right-hand side of (38). From [32, Lemma 3.2], Mignot proves the following description of the tangent cone of K 0 : withȳ a quasi-continuous representative of the function y. This provides the characterisation stated in the lemma.
The set K 0 is said to be polyhedric at (y, λ) } is the known as the polar cone. The set is polyhedric at y if it is polyhedric at (y, λ) for all λ, and it is polyhedric if it is polyhedric at y for all y. The concept of polyhedricity is useful because it is a sufficient condition guaranteeing the directional differentiability of the metric projection associated to that set (see [22,32,13] and also [42]) and this fact ultimately enables one to obtain directional differentiability for solution mappings of variational inequalities. Now, the set K 0 is not polyhedric as a subset of the space W (0, T ) since W (0, T ) lacks certain smoothness properties due to the low regularity of the time derivative. However,K 0 : The setK 0 is polyhedric as a subset of W s (0, T ) and for (y, λ) whereȳ is a quasi-continuous representative of y and which shows that (·) + : W s (0, T ) → W s (0, T ) is a bounded map. It follows that W s (0, T ) is a vector lattice in the sense of Definition 4.6 of [42] when associated to the coneK 0 . The boundedness of (·) + : W s (0, T ) → W s (0, T ) and Lemma 4.8 and Theorem 4.18 of [42] imply thatK 0 is polyhedric in W s (0, T ) and hence

Directional differentiability for VIs
We now specialise to the case where the pivot space H is a Lebesgue space, a restriction which is needed for the results of [15]. Theorem 4.1 of [15] states that for w 0 ∈ V + , the mapS w0 : where Using this, we shall first work to deduce a differentiability formula for the map S under perturbations of the right-hand side source with a fixed obstacle.
Remark 4.4. Our notation emphasises the fact that the higher-order term in (39) Let us see whyô above is a higher-order term.  (17) and (18) hold. Then belongs to L 2 (0, T ; V ) ∩ L ∞ (0, T ; H), and h is a higher-order term in L p (0, T ; H) whose convergence is uniform in d on compact subsets of L 2 (0, T ; H). The directional derivative α := ∂S z0 (f, ψ)(d) satisfies If additionally, for s ≥ 0, , (18), (20), Proof. The assumptions imply that Proposition 3.4 applies and we obtain existence of the left-hand side and the first term on the right-hand side of (43). We can via Proposition 3.11 utilise (35) with the source terms f and f + sd to write S z0 (f, ψ) and S z0 (f + sd, ψ) in terms ofS w0 , and then with the aid of the expansion formula (39), we find

Differentiability with respect to the obstacle and the source term
We clearly need some differentiability for the obstacle mapping to proceed the study further and this comes in the following assumption which we take to stand for the rest of the paper. Thus, Φ needs to be defined not just on L 2 (0, T ; V ) but on the larger space L 2 (0, T ; H). This is necessary because it implies the uniform convergence with respect to compact subsets of the direction in L 2 (0, T ; H) of the difference quotients to the directional derivative of Φ, which is a fact that we will use later in §5.5 in the analysis of some higherorder terms that arise. We could have asked for Φ to be defined on L 2 (0, T ; H) right from the outset in §1 but note that this enlargement of the domain is a restriction (and would be an unnecessary restriction for earlier sections): maps defined on L 2 (0, T ; H) are also defined on L 2 (0, T ; V ) but the converse is not true.
Using this and the expansion formula (41) forS, the second term on the right-hand side of (56) becomes where the second equality holds since every term insideS w0 on the left-hand side is in L 2 (0, T ; H) and so (39) applies. Now, plugging (57) and (58) into (56) we find where for the final equality, we again used the formula (35) which is applicable because (52) implies that z 0 ≤ Φ(ψ)(0), and we used the relation (44) between the directional derivatives ofS and S: We then set α := Φ ′ (ψ)(ρ) + ∂S(f, ψ)(d − LΦ ′ (ψ)(ρ)). From (40) where w =S(LΦ(ψ)−f ) = Φ(ψ)−S(f, ψ). Recalling the definition of α and making the substitution ϕ := Φ ′ (ψ)(ρ)+z in the above variational formulation for δ yields the formulation for α stated in the proposition. If additionally (D s ), (54), (53), then Proposition 3.11 gives the stated regularity. Dropping now the dependence on the base points for clarity, we estimate the remainder term r (which is defined as in the statement of this proposition) as follows, making use of (42) and Dividing by s and taking the limit s → 0 + , we see that the remainder term vanishes in the limit due to assumption (55). Furthermore the convergence to zero is uniform in d on compact subsets since d appears only in the final term which we know has the same property as S(·, ψ) is Hadamard differentiable. The assumption h(0) = 0 in the proposition implies that all assumptions that hold for the perturbed data also hold for the non-perturbed data (i.e. at s = 0). Without this assumption, we would have to assume in addition Φ(ψ) ∈ W r (0, T ) and (20) and (D1) along with (53).

Directional differentiability
Fix f, d ∈ L 2 (0, T ; H). We begin by choosing an element of P(f ) with sufficient regularity as described in the following assumption.
Picking u ∈ P u0 (f ) satisfying Assumption 5.1, define the sequence u n s := S u0 (f + sd, u n−1 s ) for n = 1, 2, 3, ..., u 0 s := u. Our aim will be to apply Theorem 3.8 or Theorem 3.10 to this sequence in order to show that, under additional assumptions, it is well defined and has the right convergence properties. Furthermore, we also want to obtain expansion formulae for each u n s .

Expansion formula for the VI iterates
The following sets of assumptions are to ensure that Theorem 3.8 (or Theorem 3.10) can be applied for our sequence {u n s } n∈N .
if ρ ∈ L 2 (0, T ; V ), and h : (0, 1) → L 2 (0, T ; H) is a higher-order term, then as s → 0 + , and either , v ≤ 0 for all v ∈ V + and for a.e. t, To ease the notation on the higher-order terms, we did not write the base point u in the r term (which originates from Proposition 4.9) above. We now give two results (with varying assumptions) in the next proposition concerning convergence behaviour and an expansion formula for the sequence {u n s }.
where (1) α n satisfies the VI Proof. Let us first show monotonicity of the sequence assuming existence. First take d ≥ 0. Then since u is a subsolution, u ≤ S(f, u) ≤ S(f + sd, u) = u 1 s . By the comparison principle, we find again that u 1 s and in this we obtain that {u n s } is an increasing sequence. Likewise if d ≤ 0, the sequence is decreasing. 1. FIRST CASE. Observe that since u ≤ Φ(u) and Φ(u) ∈ W s (0, T ), by (O3a), we can take the trace at t = 0 to get Φ(u)(0) ≥ u 0 . Since (D s ), (L1a) and (60) hold and as (20) is satisfied for the obstacle u (thanks to (L1a) and (O3a)), Proposition 3.4 implies that u 1 s = S(f + sd, u) ∈ W s (0, T ) exists and is increasing in time and non-negative. Regarding the upper bound for the initial data in terms of u 1 s , we argue as follows. For d ≥ 0, we may apply Φ to the inequality u 1 s ≥ S(f, u) = u and use the regularity offered by (O3a) to obtain u 0 ≤ Φ(u)(0) ≤ Φ(u 1 s )(0). If d ≤ 0, we use the condition (61) to obtain the same conclusion. Then making use of (O1a), (O3a), and (O4a) we apply Proposition 3.4 to obtain the existence for each u 2 s and subsequently, using these arguments, existence for each u n s . We now show the expansion formula (65) by induction.
Regarding the expression for the derivative, we know that ) solves the VI that appears in Proposition 4.5, i.e., as desired.
2. SECOND CASE. We will not repeat some of the same techniques used in the above case and simply focus on the differences under the different set of assumptions. Due to (L1b), u 1 s exists by Proposition 3.4. Using (62) we find u 0 ≤ Φ(u 1 s )(0). The monotonicity of {u n s } n∈N and (62) shows this bound for all u n s . Using (O2b), (O4b) and (62), we infer the existence for all u s n by the same proposition. We prove the remaining claims again by induction. For the base case, we can expand u 1 s = S(f + sd, u) by using (L1b) and Proposition 4.5 directly and we obtain a δ 1 := ∂S(f, u)(d) ∈ L 2 (0, T ; V ) ∩ L ∞ (0, T ; H) such that (67) holds. Furthermore, u 1 s ∈ W r (0, T ). Now assume the statement is true for n = k so that (68) holds. Since u k s ∈ W r (0, T ), LΦ(u k s ) ∈ L 2 (0, T ; H), and by (L2), since we have LΦ ′ (u)(α k ) ∈ L 2 (0, T ; H) and so assumption (51) of Proposition 4.9 holds, and as does (52) as shown above, and the proposition can applied to give Note that we used the fact that (O4b) implies the increasing property of all obstacles considered in the proof.
3. CONCLUSION. The claim of the VI satisfied by the α n follows from Proposition 4.9 whilst the convergence behaviour stated in the result is a consequence of either Theorem 3.8 (if d ≥ 0) or Theorem 3.10 (if d ≤ 0), using the fact that (59) implies that u 0 s is either a subsolution or supersolution.
Remark 5.5. Everything up to the convergence of the {u n s } stated in the above result holds if we do not assume (59) and either (O2a) or (O3b) respectively. Also, (L3) was necessary only to prove that each o n is a higher-order term.

Properties of the iterates
In this section, we give some basic attributes of the directional derivatives α n and the higher-order terms o n . One should not forget that these objects are time-dependent, and we will always denote the time component by t; this should not be confused with the perturbation parameter s. Lemma 5.6. The following properties hold: q.e. on {u(t) = Φ(u)(t)} for n > 1. 4. Φ ′ (u)(α n ) has the same sign as d.

The sequences
Proof. The first claim follows from the set that the α n belong to and the characterisation of Lemma 4.1. The second claim is true since the sequence u n s is increasing or decreasing in n and due to (65) and the vanishing behaviour of s −1 o n (s) and the fact that u n s − u has the same sign as d (since if d ≥ 0, u n s is increasing and hence greater than u whilst if d ≤ 0, u n s is decreasing and smaller than u). For the third claim, in (65), take the trace t = 0 (which is valid since u n s , u were defined to have trace u 0 at t = 0) to obtain Finally, we have that where the limit is in L p (0, T ; H), and hence, passing to a subsequence, the limit also holds at almost every time strongly in H: which is either greater than or less than zero depending on the sign of α n which in turn depends on the sign of d (see part 2 of this lemma).
The first part of the previous result tells us about the quasi-everywhere behaviour of the directional derivatives on the coincidence set. We can say a little more about them in an almost everywhere sense. If Φ is a superposition operator, then for each n, On the set {u = Φ(u)}, we get sα 1 ≤ −o 1 (s) and dividing here by s and sending to zero, we see by the sandwich theorem that if d ≥ 0, α 1 = 0 on {u = Φ(u)}. If Φ is a superposition operator, observe that if t is such that u(t) = Φ(u(t)), then in fact u(t) = Φ n (u(t)) for any n ∈ N.
Using this fact on the right-hand side of (69) gives us the result.

Uniform bounds on the iterates
We give a result on the boundedness of the directional derivatives α n under two different sets of assumptions. The first set requires some boundedness conditions on the obstacle mapping including a smallness condition, whilst the second requires instead some regularity and complementarity (for the latter, see [21, §7.3.1] for the parabolic VI case) for the system.
where all constants are independent of n.
Defining a n := α n L 2 (0,T ;V ) , this reads which we write as where we have denoted Solving this recurrence inequality leads to We evidently need A 2 < A 1 for this sequence to be bounded uniformly, that is, i.e., (L7a). Under this condition, the bound is uniform and there is a weak limit for a subsequence of {α n } n∈N . Since the α n are monotone, they have a pointwise a.e. monotone limit which must agree with α so indeed α n ⇀ α in L 2 (0, T ; V ).
2. UNDER REGULARITY ASSUMPTIONS. Now assume instead that (L4b) -(L6b) hold. We want to show that ϕ ≡ 0 is a valid test function in the VI (66) for α n . Thus we need to prove that The assumption (L5b) implies that −Φ ′ (u)(α n−1 ) ≥ 0 a.e. on {u = Φ(u)} and thus it belongs to T rad (36)) and this is obviously a subset of its closure in W . Therefore, 0 is a valid test function in (66) and testing with this we find which easily leads to the desired bound due to the assumption (L6b). As before, the monotonicity of the sequence implies the convergence for the whole sequence.

Characterisation of the limit of the directional derivatives
We now want to study the limiting objects associated to the sequences {α n } and {δ n }. First, we introduce some assumptions that will be of use here and in further sections.
and assume that if z : (0, 1) → L p (0, T ; H) satisfies z(s) → u as s → 0 + and b ∈ L 2 (0, T ; V ), then for all s ∈ (0, ǫ) with ǫ arbitrarily small, then where Regarding (L8), it may be helpful to note that Assumption 4.6 implies that Φ : L 2 (0, T ; V ) → W (0, T ) is completely continuous. As a precursor to characterising the directional derivative α, we study the limit of {δ n } n∈N in the next lemma. Since δ n = α n + Φ ′ (u)(α n−1 ), if Φ ′ (u)(·) : L 2 (0, T ; V ) → L p (0, T ; H) is bounded, we can find a subsequence of {δ n } n∈N such that δ nj ⇀ δ in L p (0, T ; H) for some δ. In fact under additional assumptions the convergence holds for the full sequence as shown below.
To characterise δ as the solution of a VI in itself, it becomes useful to define the set Proof. From (45), δ n satisfies the VI If (L8) holds, we can pass to the limit here using the continuity of Φ ′ (u)(·) : L 2 (0, T ; H) → W (0, T ) and the limiting object δ satisfies the inequality stated in the lemma. We must check that It is clear that the orthogonality condition is satisfied due to the convergence in the previous lemma. By (37),

Due to Mazur's lemma, there is a convex combination
Since this convergence is strong, for a subsequence, v km (t) → δ(t) in V and hence pointwise q.e. for a.e. t ∈ [0, T ].
By definition, δ n (t) ≥ 0 everywhere on and using the fact that a countable union of capacity zero sets has capacity zero and the inequality (70), we can pass to the limit to deduce that δ(t) ≥ 0 quasi-everywhere on {u(t) = Φ(u)(t)} for a.e. t ∈ [0, T ].
Proposition 5.13. Under the conditions of the previous lemma, α satisfies the QVI Proof. From the definition of α n in terms of δ n in (63), we obtain Using this fact in the QVI for δ given in Lemma 5.12 yields which translates into the desired result after setting w := Φ ′ (u)(α) + z.

Dealing with the higher-order term
We come to the final part which consists in showing that the limit of the higher-order terms o n is indeed a higher-order term itself. The idea is to be able to commute the two limits We estimate the third term above using (71) for some C < 1 by the assumption (L13). The above can be recast as The recurrence inequality can be solved for a n in terms of a 1 and the b i : a n ≤ C n−1 a 1 + C n−2 b 1 + C n−3 b 2 + ... + Cb n−2 + b n−1 .
By the weak convergence of the α n in L 2 (0, T ; V ) (and hence strong convergence in L 2 (0, T ; H)), {α n−1 } n∈N is a compact set in L 2 (0, T ; H). Since Φ : L 2 (0, T ; H) → W (0, T ) is Hadamard differentiable, it is compactly differentiable, meaning that s −1 l(s, α n ) → 0 in W (0, T ) uniformly, thus the first three terms in the definition of b n are taken care of. For the final term, we note that {LΦ ′ (u)(α n−1 )} n∈N is also compact in L 2 (0, T ; H) by (L9). Therefore, for any ǫ > 0, there exists a τ 1 > 0 independent of j such that STEP 4. As o 1 (s) = r(s, 0, 0) = o(s, d) is a higher-order term, we know that there is a τ 2 > 0 such that Now recalling (73), for s ≤ min(τ 1 , τ 2 ), (74) and (75) This shows that o n (s)/s tends to zero uniformly in n.
We are finally in position to state the main theorem which condenses the various intermediary results we obtained above.
Let us relate this result to notions of differentiability commonly used in set-valued and multi-valued analysis. Recall that the map P : F ⇒ U has a contingent derivative β at (f, u) (with u ∈ P(f )) in the direction d, written β ∈ DP(f, u)(d), if there exist sequences β n → β, d n → d and s n → 0 such that u + s n β n ∈ P(f + s n d n ).
We claim that P does indeed possess contingent derivatives and our main results furnishes us with one such contingent derivative. Proof. Indeed, given u ∈ P(f ) and sequences s n → 0 and d n ≡ d, we can, by Theorem 5.15, find u s n ∈ P(f + s n d) such that u s n = u + s n α n + o(s n ; d) and hence u + s n β n = u s n ∈ P(f + s n d n ), where the sequence β n := α n + o(s n ; d) s n = u s n − u s n is such that β n → α as n → ∞ because the term o(s n ; d) is higher order.

Other approaches to differentiability
We now discuss some possible alternative approaches to deriving Theorem 5.15 (or a variant of this theorem). The idea here is to bootstrap by using the already-achieved results on elliptic QVIs in [2], however, we shall see that this is not at all straightforward.

Time discretisation and elliptic QVI theory
One idea is to apply the elliptic QVI theory of [2] to the time-discrete problems and then pass to the limit there. With the definition Ψ n,N (z) := Φ(z)(t N n ), recall the notations defined in §2; in particular Q t N n is the solution mapping defined as z ∈ Q t N n (g) if and only if z ∈ V, z ≤ Ψ n−1,N (z) : With u N 0 := u 0 ∈ V for a given initial data, set Under the assumptions on the source f and the direction d in §2, we know by [2,Theorem 1] In a similar way, since is a higher-order term and since u N Along these lines, we obtain the following lemma.
Lemma 6.1. Given u N n defined as above, there exists u N s,n ∈ Q t N n−1 (hf N n + u N n−1 + shd N n ) and γ N n such that and o N n (s) = o N n (s, hd N n + γ N n−1 , o N n−1 (s); hf N n + u N n−1 ) is a higher-order term. Here, A(u N n ) = {}. Proof. We prove this by induction. The base case has been shown above. Suppose it holds for the nth case. Then given (thanks to the formula relating u N s,n and u N n ) such that (79) It is then clear from (79) (since the structure of the time-discretised QVI is the same as those considered in §2 and we can apply the same argumentation as there) that the interpolants u N s ,û N s constructed from {u N s,n } n∈N and which satisfy are bounded in the appropriate spaces and hence there exists a u s such that u N s → u s and i.e., u s ∈ P(f + sd). Clearly, similar claims are true for u N ,û N constructed as interpolants from {u N n } n∈N . We now need to obtain uniform bounds on the γ N n and the higher-order term. Lemma 6.2. Suppose that v = 0 is a valid test function in (78) (which is the case in the VI setting). Then and hence γ N L ∞ (0,T ;H) + γ N L 2 (0,T ;V ) ≤ C. Proof. Since 0 is feasible, the proof for the boundedness of the first two estimates follows that of the first two estimates of Lemma 2.3. The proof of Corollary 2.4 shows that these bounds imply the boundedness in the Bochner spaces above.
Multiplying (77) by χ T N n−1 and summing up over N , we find u N s = u N + sγ N + o N (s). Due to the bounds on the left-hand side and the first two terms on the right-hand side, we can pass to the limit in o N too and we find u s = u + sγ * + o * (s), for some γ * ∈ L 2 (0, T ; V ) ∩ L ∞ (0, T ; H), and the equality also holds in this space. It remains to be seen that o * is a higher-order term and to characterise the term γ * .
Here we come to a roadblock since we do not know how the terms o N depend on the moving base point and thus we cannot say anything about the uniform convergence of the o N and are unable to identify the limit of the higher-order terms as higher order; recall that we warned the reader of this issue in Remark 4.4. Alternatively, if we were able to identify γ * as α from the previous sections, this would suffice to show the higher-order behaviour. But it seems difficult to handle the constraint set satisfied by γ N n and obtain a type of Mosco convergence for recovery sequences to approximate the limiting test function space.

Elliptic regularisation of the parabolic (Q)VI
Another possible approach is through regularising the parabolic QVI as an elliptic QVI, applying the elliptic theory of [2] to the regularised problem and then passing to the limit in the regularisation parameter. This elliptic regularisation of parabolic problems can be seen in the work of Lions [30, p. 407].
The idea is to include in the parabolic inequality a term involving the second time derivative of the solution, i.e., −ǫu ′′ with ǫ > 0; more precisely we wish to consider the QVI (80) The 'final time' condition is necessary to have a well defined problem. Define V := {v ∈ W s (0, T ) : v(0) = 0}. Integrating by parts in (80) and using the initial and final conditions on u and the test function space, we obtain the weak form: find u ∈ V such that (81) Thinking of the time component as simply another space dimension, the resulting elliptic operatorÂ : which is clearly T-monotone, bounded and coercive in the Sobolev-Bochner space V; for the latter, observe that Clearly the solution of (81) depends on ǫ so let us write it as u ǫ . Working formally, we may apply [2, Theorem 1.6] and we find existence of a u ǫ s solving (81) with source term f + sd and a α ǫ such that where o ǫ is a higher-order term. As for α ǫ , it satisfies Testing (81) with v = 0, we obtain the bound In a similar way, u ǫ s is also bounded in this space uniformly in ǫ and s. Let us also suppose that we have a uniform bound on α ǫ . Then it remains to be seen that the limit of the o ǫ is a higher-order term, and this is where we once again run into problems. A monotonicity in ǫ result on the solutions of (81) would be useful.
It is worth remarking that this technical issues also arises in the VI case and in the simpler unconstrained (PDE) case.

Example
We study here an example where the obstacle mapping is given by the inverse of a parabolic differential operator. This example is motivated by applications in thermoforming (which is a process that manufactures shapes such as car panels from a given mould shape; see [2] and references therein for more information): in [2], we studied an elliptic nonlinear PDE-QVI model that describes such a thermoforming process but in reality, the full model is a highly complicated evolutionary system of equations and QVIs involving obstacle maps that are inverses of differential operators. As a first initial step on the road to studying the full problem, we consider a simplification and study the following example.
We denote by C B b and C B a the constants of boundedness and coercivity for the operator B while A : V → V * denotes an operator which is assumed to satisfy all assumptions given in the introduction of this paper and it should also give rise to a bilinear form which is smooth. Letting W = H 2 (Ω), we further assume that This equivalence of norms assumption is related to (boundary) regularity results which follow given enough smoothness of the boundary and depending on the specific elliptic operators, see [16,Theorem 4,§6.3]. The equivalence of norms implies Define Φ(ψ) := w as the solution of (82). By the standard theory of parabolic PDEs, Φ : L 2 (0, T ; H) → W s (0, T ). This regularity implies (using the equation itself) that for ψ ∈ L 2 (0, T ; H), BΦ(ψ) ∈ L 2 (0, T ; H) and hence by (83) and the assumption on the range of A, Φ(ψ) ∈ L 2 (0, T ; W ) and AΦ(ψ) ∈ L 2 (0, T ; H). It follows then that LΦ(ψ) ∈ L 2 (0, T ; H) too (recall that L = ∂ t + A). Since Φ(ψ) ∈ L 2 (0, T ; W ) ∩ H 1 (0, T ; H), the theorem of Aubin-Lions [14,Theorem II.5.16] gives the continuity in time Φ(ψ) ∈ C 0 ([0, T ]; V ).
We proceed now with checking the assumptions that will eventually enable us to use Theorem 5.15.
Lemma 7.1. Assumption 4.6 on the Hadamard differentiability of Φ holds and the directional derivative satisfies where Φ 0 is the solution mapping associated to the equation (82) with zero initial condition (i.e. Φ 0 is the same as Φ except for the initial condition).
Given a source term f ∈ L 2 (0, T ; H + ) and initial condition u 0 ∈ V + satisfying u 0 ≤ Φ(0)(0) = w 0 and (60), we fix a solution u solving the QVI Find u ∈ W (0, T ) with u ≤ Φ(u) : (84) The fact that such a solution indeed exists will be shown later. Proof. We showed in the previous paragraph the satisfaction of (O3a), and the condition (61) is irrelevant since d ≥ 0. Let us check the remaining assumptions.
• Regarding the increasing property of Φ in time, we argue as follows. Making the substitutionw = w − w 0 in (82) in order to remove the initial condition, the transformed PDE is (g(ψ(t 1 − r, y)) − g(ψ(t 2 − r, y))G(x, y, r) dydr − t2 t1 L 0 (g(ψ(t 2 − r, y)) − Bw 0 (y))G(x, y, r) dydr ≤ 0 since ψ is increasing in time and g is increasing which implies that the first term is less than zero, and the second term is also clearly non-positive as its integrand is non-negative (due to the assumption on w 0 ). This shows that t →w(t) = w(t) − w 0 = Φ(ψ)(t) − w 0 is increasing, implying (O4a).
• Multiplying the equation by w − leads to since g ≥ 0, which shows that Φ(ψ) is non-negative too as long as w 0 ∈ V + , thus (O1a) and (L1a).
Subtracting one from the other and testing with the difference, we obtain Making use of Young's inequality on the right-hand side and then integrating over time, we find due to the continuity of g in L 2 (0, T ; H) that w n → w in L 2 (0, T ; V ) ∩ L ∞ (0, T ; H) and so (O2a) is also valid.
• From the standard energy estimate i.e., (L4a) after using the continuity of the embedding V ֒→ H.
• The second estimate above also leads to • For (L5a), we simply use the equation itself and (85): Therefore, if T and/or g ′ ∞ is sufficiently small, assumption (L7a) holds. • The estimate (85) yields the first property in (L3) whilst (85) and (87) yield the second property.
• We also have that LΦ ′ (u)(h n ) = (∂ t + A)Φ 0 (g ′ (u)h n ) = (∂ t + A)w n = (∂ t + B)w n + Aw n − Bw n , which implies that The first term on the right-hand side converges to zero because g ′ is bounded. For the third term, we note the following standard estimate for linear parabolic PDEs (using the differentiability of B): Regarding the second term of (88), we manipulate using the equivalence of norms as following: A(w n − w) L 2 (0,T ;H) ≤ C 1 w n − w L 2 (0,T ;W ) ≤ C 2 w n − w L 2 (0,T ;H) + B(w n − w) L 2 (0,T ;H) and we apply again the estimate from the previous step on the right-hand side and this eventually leads to (L9).
Having met all the requirements, we may apply Theorem 5.15 to infer that the solution mapping taking f → u in (84) is directionally differentiable in the stated sense and the associated derivative solves the QVI