On an Anisotropic Fractional Stefan-Type Problem with Dirichlet Boundary Conditions

In this work, we consider the fractional Stefan-type problem in a Lipschitz bounded domain $\Omega\subset\mathbb{R}^d$ with time-dependent Dirichlet boundary condition for the temperature $\vartheta=\vartheta(x,t)$, $\vartheta=g$ on $\Omega^c\times]0,T[$, and initial condition $\eta_0$ for the enthalpy $\eta=\eta(x,t)$, given in $\Omega\times]0,T[$ by \[\frac{\partial \eta}{\partial t} +\mathcal{L}_A^s \vartheta= f\quad\text{ with }\eta\in \beta(\vartheta),\] where $\mathcal{L}_A^s$ is an anisotropic fractional operator defined in the distributional sense by \[\langle\mathcal{L}_A^su,v\rangle=\int_{\mathbb{R}^d}AD^su\cdot D^sv\,dx,\] $\beta$ is a maximal monotone graph, $A(x)$ is a symmetric, strictly elliptic and uniformly bounded matrix, and $D^s$ is the distributional Riesz fractional gradient for $0<s<1$. We show the existence of a unique weak solution with its corresponding weak regularity. We also consider the convergence as $s\nearrow 1$ towards the classical local problem, the asymptotic behaviour as $t\to\infty$, and the convergence of the two-phase Stefan-type problem to the one-phase Stefan-type problem by varying the maximal monotone graph $\beta$.


Introduction
The classical Stefan problem, in an open bounded Lipschitz domain Ω x = (x 1 , . . . , x d ) and for time t ∈ [0, T ], can be formulated in Q T = Ω×]0, T [ by an evolution equation involving a subdifferential operator where ϑ(x, t) is the temperature, D is the gradient, A = A(x) is a symmetric, strictly elliptic and bounded matrix, and β corresponds to a maximal monotone graph, such that β(r) = b(r) + λχ for χ ∈ H(ϑ) for the maximal monotone graph H(r) associated with the Heaviside function, i.e. H(r) = 0 for r < 0, H(r) = 1 for r > 0, H(0) = [0, 1], and b a given continuous and strictly increasing function, λ > 0 (see Figure 1) with inverse γ = β −1 satisfying lim r→+∞ γ(r) = +∞ and lim r→−∞ γ(r) = −∞ for the two-phase problem and γ(r) = 0 for r ≤ λ for the one-phase problem. The notation β(ϑ) should be understood as follows: there exists a section η of the multifunction β(ϑ) which satisfies the required conditions. In turn, ϑ is easy to recover from η since β −1 = γ is a single-valued mapping. For works on the variational formulation of the classical Stefan problem, see for instance [38], [31], [28], Chapter V.9 of [33], Section 3.3 of [35], [19], [50], [40], [41] and [52]. We can also consider the one-phase problem (I) as the limit of the two-phase problem (II). Indeed, physically, for large Stefan number, the liquid phase only contributes exponentially small terms to the location of the solid-melt interface. Therefore, at times close to complete solidification, the temperature in the liquid essentially vanishes and the two-phase problem reduces to the one-phase problem. For more detailed discussions, see [37]. See also [49] for the one-dimensional case in the classical setting s = 1.
Here, we consider the corresponding fractional Stefan-type problem, given in Q T by where L s A = −D s · AD s is a non-local operator defined with the distributional Riesz fractional derivatives, with anisotropy given by a measurable matrix A = A(x), which is symmetric, strictly uniformly elliptic and bounded independent of time satisfying a * |z| 2 ≤ A(x)z · z ≤ a * |z| 2 (1. 3) for almost every x ∈ R d and all z ∈ R d . Then, the classical problem (1.1) corresponds to the case s = 1, i.e. (1.2) with the operator L 1 A , where D 1 = D. The operator L s A can be viewed as an anisotopic generalisation of the fractional Laplacian. Indeed, following the works of Silhavy [47], Shieh-Spector [45]- [46] and Comi-Stefani [16]- [17], the Riesz fractional s-gradient (D s ) and the s-divergence (D s ·) are defined in integral form for sufficiently regular functions u and vector fields φ, respectively, by . Then, in the distributional sense, it is well-known that −D s · D s u = (−∆) s u for u ∈ C ∞ c (Ω) (see for instance, [45], [47]), where (−∆) s is the fractional Laplacian defined as |x − y| d+2s dy for 0 < s < 1.
Furthermore, we have the convergence of the fractional derivatives to the classical derivatives as s 1, i.e. D s u → Du, as in Comi-Stefani [16], Bellido et al. [7] and Lo-Rodrigues [36], for u ∈ H 1 (R d ) and u ∈ H 1 0 (Ω). In this work, we are concerned with the classical fractional Sobolev space H s 0 (Ω) in a bounded domain Ω ⊂ R d with Lipschitz boundary, for 0 < s < 1, defined as where u is extended by 0 in R d \Ω, so that this extension is also in H s (R d ). By the classical fractional Poincaré inequality (see Lemma 2 below), we shall consider the space H s 0 (Ω) with the following equivalent norm u for 1 ≤ q < 2 * , where 2 * = 2d d−2s and q > 2 # = 2d d+2s when s < d 2 , and if d = 1, 2 * = q for any finite q and 2 # = q q−1 when s = 1 2 and 2 * = ∞ and 2 # = 1 when s > 1 2 . We recall that those embeddings are continuous also for q = 2 * when s < d 2 (see for example, Theorem 4.54 of [25]). The nonlocal operator L s A = −D s · AD s may be defined in the duality sense for u ∈ H s (R d ): L s A u, v :=ˆR d AD s u · D s v ∀v ∈ H s 0 (Ω), (1.6) with v extended by zero outside Ω, defining an operator from H s (R d ) to H −s (Ω) since AD s u ∈ L 2 (R d ) d . Also for u ∈ H s 0 (Ω), since we can extend it by 0 outside Ω to obtain a function in H s (R d ), L s A : H s 0 (Ω) → H −s (Ω) can also be represented by L s A u = −D s · (AD s u). (1.7) Given anyg ∈ H s (R d ), we introduce g ∈ H s (R d ) defined on the whole space R d which satisfies g| Ω c =g and is L s A -harmonic in Ω, that is to say, we solve the Dirichlet problem with g =g a.e. on Ω c for the equation in a weak sense, which meansˆR d AD s g · D s v = 0 ∀v ∈ H s 0 (Ω).
Remark 1. Note that the variational problem (1.12) incorporates the Dirichlet condition (1.11) in the original problem given in (1.2) because of the definition (1.8). Since this implies´R d ×[0,T ] AD s g · D s ξ = 0 for all ξ ∈ Ξ s T , we obtain (1.14) without that term. Although in general η, ϑ may be nonzero outside Ω, except for the bilinear form´R d ×[0,T ] AD s ϑ · D s ξ, the other integral terms in the variational formulation (1.14) are only integrated over Ω in space, since the test function ξ is 0 in Ω c ×]0, T [. Different non-local versions of Stefan-type problems have previously been considered, including in [9] and [15] for nonsingular integral kernels, in [53], [8], [44] and [42] for the fractional Caputo derivatives, and in [21], [22], [23], [24] and [30] for the fractional Laplacian and its nonlocal integral generalization in [2]. Stefan-type problems that are fractional in the time derivative have also been considered (see, for instance, [43], [34] and [14].) Indeed, when the matrix A is a multiple of the identity matrix, the fractional Stefan-type problem (1.2) reduces to that with the fractional Laplacian as considered in [21]- [24]. Furthermore, in instances as described in Section 2.3 of [36] when the fractional operator L s A is replaced with a nonlocal operatorL s a , corresponding to a Dirichlet form with the kernel a which satisfies some compatibility conditions, (1.2) may also be considered a nonlocal Stefan problem, as considered in [2]. However, an equivalence relation between the fractional operator with the matrix A and the nonlocal operator with the kernel a cannot be established in general except in the isotropic homogeneous case (for more details, see Section 2.3 of [36]), so the two Stefan-type problems with those two operators are not equivalent.
In this paper, we show the existence of a unique solution for the fractional Stefan-type problem with Dirichlet boundary conditions, where the spatial operator is a general anisotropic non-local singular operator of fractional type as given by (1.6), and we keep the classical temperature-enthalpy relation illustrated in Figure 1. This relation in the classical equation (1.1) incorporates, in a generalised form, the free boundary condition relating the balance between the normal velocity of the interface and the jump of the local anisotropic heat flow. In 1-dimension, the extension of the classical free boundary Stefan condition to fractional diffusion, as in the recent paper [44] with the fractional Caputo derivative in the nonlocal diffusive term, can be easily made explicit. Similar explicit formulation can be used with the 1-dimensional fractional Riesz spatial derivative when, for each fixed time, the free boundary is a point.
However, in higher dimensions, the Riesz fractional s-gradient, as proposed in [47], is an appropriate fractional operator maintaining translational and rotational invariance, as well as homogeneity of degree s under isotropic scaling, and so the L s A operator gives a natural and appropriate anisotropic generalisation of the fractional Laplacian. Keeping the generalised Stefan condition in the evolution equation (1.2) involving the maximal monotone operator β is a natural generalisation for the formulation of the anisotropic Stefan problem, extending [23] and [24], which corresponds to the case where the matrix A is the identity matrix in the unbounded domain. Such an anisotropic operator is coordinate invariant, which makes it more suitable in higher dimensions. Furthermore, the use of this L s A operator allows us to recover the classical Stefan problem when s = 1, which is in accordance with a requirement of weak continuity from the nonlocal model to the local model, when s 1. However, a main issue remains open in the fractional multidimensional model, namely what is the physical meaning of the Stefan condition due to the lack of a convenient interpretation and definition for the fractional heat flux across the solid-liquid interface.
In Sections 2 and 3, we employ Hilbertian techniques to show the existence of a generalised enthalpy solution and a weak temperature solution to the initial and boundary value two-phase Stefan-type problem (1.12)-(1.14) following the approach of Damlamian [18]- [19] for the classical case s = 1.
Making use of convergence properties of the fractional derivatives to the classical derivatives when s 1, we show, in Section 4, that the solution of the fractional Stefan-type problem converges to the solution of the classical case corresponding to s = 1. Next, we consider the asymptotic behaviour of the solution as t → ∞ in Section 5. Such convergence properties apply to both the two-phase problem, and the onephase problem, which corresponds to the case of a nonnegative temperature. The one-phase problem (I) is recovered in Section 6 from the two-phase problem (II), when the maximal monotone graph for (II) (see Figure 1) degenerates to that of the one-phase problem (I).
Finally, we complete our paper with three appendices: A on the time dependent Dirichlet problem for the fractional operator L s A , B on the variational inequality formulation for the two-phase and the one-phase problems, and C on the stability of the eigenvalues and eigenfunctions for the operator L s A in H s 0 (Ω) with respect to s including the convergence s 1.

Existence of the Generalised Enthalpy Solution η
Let L = L s A be the duality mapping defined by from H s 0 (Ω) to H −s (Ω) with H s 0 (Ω) identified to a subspace of L 2 (Ω). Here ·, · is the duality between H −s (Ω) with H s 0 (Ω), with u, v extended by zero outside Ω. The equality of the inner product in H −s (Ω) given by (·, ·), with the topology endowed from L, with the equivalent inner product [·, ·] A in H s 0 (Ω) holds by Riesz representation theorem, with U = Lu and V = Lv respectively. (2.2) This is possible by assumption (1.3) and the Poincaré inequality, as long as Ω is bounded.
In this section, we consider the two-phase problem with γ is Lipschitz with Lipschitz constant C γ such that γ(0) = 0 and lim inf |r|→+∞ γ(r) r > 0.
We prove an existence theorem for the enthalpy η similar to the classical case, as given in [19] and [18] (See also [52] for further developments). To do so, we need a result of Attouch-Damlamian [5]- [6] in the case where the Hilbert space H is H −s (Ω).

Proposition 1.
[Theorem 1 of [5], and [6]] Let (ϕ t ) t∈[0,T ] be a family of lower semi-continuous convex functions on a Hilbert space H. Assume that there exists a function a belonging to BV (0, T ) such that the following holds: Furthermore, the following estimates hold independent of ϕ: Making use of this proposition, we can show the following existence result.
Proof. We apply Proposition 1 with |a(t) − a(τ )| = g(t) − g(τ ) L 2 (Ω) to the following functions φ t on the Hilbert space H −s (Ω) given for each t ∈ [0, T ] by where j is the primitive of γ such that j(0) = 0. Then, j is quadratic and the domain D(φ t ) of φ t is given by thanks to the Cauchy-Schwarz inequality and making use of the fact that W lies in L 2 (Ω). It is well-known (see for instance, Theorem 17 of [11]) that φ t is lower semi-continuous, convex, proper and coercive on H −s (Ω). Furthermore, there exist constants δ and c such that Consequently, by classical results of subdifferentials (see for instance, [11] or [32]), the subdifferential ∂φ t is a maximal monotone operator of H −s (Ω). In fact, the subdifferential ∂φ t is characterized as follows: if and only if U ∈ L 2 (Ω) and L −1 (V ) + g = γ(U ) a.e. in Ω, (2.17) and we recall from (2.14) that , representing the Dirichlet condition in weak form in the trace sense for s > 1 2 and more generally γ(U ) = g in Ω c . Indeed, the characterisation of the subdifferential in terms of the convex conjugate functions involving (U, V ) for U, V ∈ H −s (Ω) reads as: Set J(W ) =´Ω j(W ). Recognising the evaluation at the point L −1 V + g with the convex conjugate on L 2 (Ω) of j(W ), by well-known results (see for example Lemma 1 of [12], or [39]), we can associate the convex conjugate J * (U ) with´Ω j * (U ), so we have where j * is the convex conjugate of j on R. From (2.18), this means that Recall (see for example, [4]) that for dual convex functions j and j * , for all numbers a, b. Therefore, the integrand in (2.19) must be non-negative, and so it is almost everywhere zero, i.e.
This means that the points U and L −1 V + g are conjugated, i.e. L −1 V + g ∈ ∂j(U ). By definition of j as the primitive of γ, we have L −1 V + g = ∂j(U ) = γ(U ). Now, we are ready to apply Proposition 1 in the space H −s (Ω) with the convex functions φ t . For W ∈ D(φ τ ) ∩ D(φ t ) = D(φ 0 ) since the domain D(φ t ) as given in (2.15) is independent of t, we have, by (2.14), so, by the Cauchy-Schwarz inequality, (2.20) Also, from (2.16), we have that where C 5 depends only on δ, |Ω| and g BV (0,T ;L 2 (Ω)) . Therefore, with the given regularity of g inherited fromg, (2.4) is satisfied, hence we can apply Proposition 1 to solve the Cauchy problem Moreover, the estimates in Proposition 1 and (2.21) give which is (1.12). Finally, for ξ ∈ Ξ s T , we can integrate in time by parts and obtain (1.14). Remark 2. To apply Proposition 1, we see from (2.20) that it is sufficient to require g ∈ BV (0, T ; L 2 (Ω)), as in [19]. However, we require additionally that g ∈ L 2 (0, T ; H s (R d )) so that (1.12)-(1.14) is well-defined.
Remark 3. We observe that the general result of the above proposition and theorem applies to general maximal monotone operators of subdifferential type with different functions γ, and so, besides two-phase Stefan-type problems, it applies also to other models including the porous medium equation. In fact, different assumptions on γ can be used (see page 12 of [18] for more details), generalising the case of the assumption (2.3).

Remark 4.
Considering the above proposition in the case where the Hilbert space H is H −s (Ω), the solution to the Cauchy problem (2.5) in H −s (Ω) with the convex function φ t , with domain L 2 (Ω), is obtained by considering the approximated problem with the convex function given by its Yosida approximation φ t,λ (V ) = 1 λ (Id + (Id + λ∂φ t ) −1 )V . Since the estimate (2.4) carries over to φ t,λ , we can apply the Gronwall's inequality to obtain the estimates (2.6) and (2.8) for the solutions to the approximated problem as in the Part 3 of the proof of Theorem 1 of [5]. Next, we make use of the absolute continuity of the map t → φ t,λ (V ) to apply to the (2.5) to obtain the estimate (2.7) from the time derivative. Passing to the limit for the approximated problems give the corresponding constants C 1 , C 2 and C 3 for the problem (2.5) in H −s (Ω).
Therefore, for σ ≤ s ≤ 1, recalling that we have the continuity of the inclusions H −σ (Ω) ⊂ H −s (Ω) ⊂ H −1 (Ω) as a consequence of Lemma 3 below, we can bound the H −s (Ω) norms with the H −σ (Ω) norms, thereby obtaining the solution to (2.5) for all s, σ ≤ s ≤ 1 with the corresponding estimates (2.6)-(2.8) for the constants C 1 , C 2 and C 3 depending only on σ and independent of s.
Since the constant C 4 is obtained from (2.8) and (2.21), similarly, we can once again consider the problem (1.12) in H −s (Ω) for each s, σ ≤ s ≤ 1, and such that the constant C 4 in (2.12) may be chosen depending only on σ and not on s.
Lemma 1. Let F ∈ L 1 (0, T ) and y ∈ L ∞ (0, T ) be non-negative functions and C > 0 a constant such that

Then we have
Remark 5. In general, for γ ≡γ, g =ĝ and an arbitrary time interval 0 ≤ t 1 < t 2 ≤ T , with a similar argument we have the fractional version of the continuous dependence property corresponding to Lemma 3.2 of [20] for the classical case s = 1: As a consequence, we immediately see that if f =f , g =ĝ and γ =γ, then for the same given data.

Regularity of the Weak Temperature Solution ϑ
If we further assume that g has two time derivatives, by the Lipschitz continuity of γ, we can achieve higher regularity of the weak temperature solution ϑ = γ(η) in (1.14). The proof makes use of the Faedo-Galerkin method, and follows closely Chapter 6 of [18], and we include it here for completeness.
Let (F n ) n∈N be an increasing set of finite dimensional subspaces of H s 0 (Ω), such that their union is dense in H s 0 (Ω), generated by the eigenvectors of the operator L −1 | L 2 (Ω) . This is possible since the inverse of L is compact in L 2 (Ω), by the compactness of the injection H s With this proposition, our approach would be to determine the subdifferental of φ t,n and show that they converge to φ t in the sense of Mosco. We recall that ϕ n M − → ϕ if for every x ∈ D(ϕ), there exists an approximating sequence of elements x n ∈ D(ϕ n ), converging strongly to x, such that lim sup n→∞ ϕ n (x n ) ≤ ϕ(x), and for any subsequence ϕ n k of ϕ n such that x k x in H, we have lim inf k→∞ ϕ n k (x k ) ≥ ϕ(x). Then applying Proposition 3 to our Faedo-Galerkin approximation, and with the additional estimates we obtain from Proposition 1, we can pass to the limit to get the additional regularity to the solution for the limit problem.
For simplicity, we drop the parameter t and consider t to be fixed in ]0, T [, and we denote φ t,n as φ n and φ t = φ. Denote i to be the compact injection of H s 0 (Ω) into L 2 (Ω) and take E n = i(F n ) by considering F n as a subspace of H s 0 (Ω). It is clear that i −1 is an isomorphism between E n and F n , with norm depending on n.
Indeed, for an eigenvector u of L| L 2 (Ω) corresponding to an eigenvalue µ, we have, by definition, . In addition, since γ satisfies the growth condition (2.3) at ±∞, its primitive j is quadratic at ±∞ (so that j(r)/|r| 2 and its inverse remain bounded as r → ±∞). Therefore, by the dominated convergence theorem, ). On the other hand, the sequence φ n is decreasing (since F n is increasing), so we conclude the Mosco convergence of φ n to φ given that φ is known to be lower semi-continuous.
Next, we want to obtain a solution of the approximate Cauchy problem for η n , making use of Proposition 1 as in the proof of Theorem 1.
Proof. Denote the inf-convolution of two convex functions by the composition operator ∇. Then by definition, we know that the convex conjugate φ * n = (φ * ∇I * F * n ) * * , where the double asterisk * * stands for the regularized l.s.c. function of ψ n = φ * ∇I * where the orthogonality is inherited from the duality between H s 0 (Ω) and H −s (Ω). Since L(F n ) = F * n , (F * n ) ⊥ is also the orthogonal of F n in H s 0 (Ω). We therefore have Since γ is globally Lipschitz, β satisfies the growth assumption (2.3) at infinity, so the function j * is quadratic at infinity and therefore z →´Ω j * (z) is continuous in L 2 (Ω). Furthermore, the function z →´Ω j * (z) is coercive in L 2 (Ω). Henceforth, we deduce that there exists z = z(v) in L 2 (Ω), not necessarily unique, such that ψ n (v) = Hence, taking a vector ξ in that basis, we have Since z(v) is the weak limit in L 2 (Ω), considering a minimising sequence of such i(z), we have the result.
Futhermore, using the coercivity of the integral of j * in L 2 (Ω) again, we see that ψ n is lower semicontinuous in H s 0 (Ω), so ψ n = φ * n . Therefore, V ∈ ∂φ n (U ) if and only if i * (U ) ∈ F * n , and there exists But since U ∈ D(φ) ∩ F * n ⊂ L 2 (Ω), we can rewrite this aŝ Ω j(U ) + j * (g + z) − U (g + z) = 0, so, as in the proof of Theorem 1, we have that the points U and g + z are conjugated by j, thus z(v) + g = ∂j(U ) = γ(U ) a.e. in Ω. The reverse is also clearly true.
Making use of these properties, we can therefore multiply (3.5) by ∂hn ∂t ∈ L 2 (0, T ; F n ) to obtain and from f n = P En f , we obtain Ω ∂η n ∂t ∂h n ∂t Now, recalling the definition ofh n , we observe that Integrating this over [0, t] for t ≤ T , we obtain, by the coercivity of A in (1.3) and integrating by parts in time,ˆt 0ˆΩ ∂η n ∂t ∂γ(η n ) ∂t (3.10) Now, we know by (2.8) that φ t,n (η n (t)) is bounded independent of n and t, so η n L ∞ (0,T ;L 2 (Ω)) is bounded independent of n (see also (2.10)). Then, by the Cea-type lemma (see, for instance, Proposition 2.5 of [1]) given by we have, by the compatibility of the initial condition giving h n (0) = P Fn (γ(η n (0))−g(0)) = P Fn (ϑ(0)−g(0)), Now, letting C γ be the Lipschitz constant of γ, we have Also, observe the boundedness ofη n in L ∞ (0, T ; L 2 (Ω)), since η n is obtained as a solution to the Faedo-Galerkin finite dimensional approximated problem (3.3) and therefore also satisfies (2.12). Therefore, applying the Cauchy-Schwarz inequality to the term´t 0´Ω f n ∂γ(ηn) ∂t and making use of the assumption Using again (3.11), we can easily take the supremum over all time to obtain the second uniform bound P Fn (γ(η n ) − g) L ∞ (0,T ;H s 0 (Ω)) = h n L ∞ (0,T ;H s 0 (Ω)) ≤ C 7 . Remark 6. For fixed σ > 0 and s such that σ ≤ s ≤ 1, similarly to Remark 4, we observe thatη n ∈ L ∞ (0, T ; L 2 (Ω)), andη n can be bounded for each s by a constant depending on σ but independent of s, by the continuity of the eigenfunctions (in Appendix C), and depending explicity on T and γ. Similarly, by Appendix C, the η n 's are bounded independent of s ≥ σ in H 1 (0, T ; F * n ). This allows us to consider the convergence of the variational problem as s varies.

Remark 7.
It can be seen that the bounds in (3.4) can be made to depend only on σ > 0 and independent of s for σ ≤ s ≤ 1, by the continuity of the eigenfunctions as shown in Appendix C. Then, as in Remark 6, the bounds D s (ϑ − g) L ∞ (0,T ;L 2 (R d ) d ) and ∂ϑ ∂t L 2 (0,T ;L 2 (Ω)) in (3.13) only depend only on σ and independent of s, allowing us to consider the convergence of the variational problem as s varies. 4 Convergence to the Classical Problem as s 1 Next, as s 1 the s-fractional derivatives converge to the classical derivatives, we show that the corresponding solutions to the fractional Stefan-type problem converge in appropriate spaces to the classical one. We first recall the fractional Poincaré inequality.
Lemma 2 (Fractional Poincaré inequality, Theorem 2.9 of [7]). Let s ∈ (0, 1). Then there exists a constant C P = C(d, Ω) > 0 such that for all u ∈ H s 0 (Ω). To consider the convergence of the problem as s 1, we start with a continuous dependence property of the Riesz derivatives as s varies, which can be easily shown using Fourier transform first for u(t) ∈ C ∞ c (Ω), and extended by density as in Lemma 3.7 of [36].
s ] for 0 < σ < s ≤ 1. As a consequence, we have the following estimate: for σ ≤ s ≤ 1, for any u(t) ∈ H s 0 (Ω) for a.e. t ∈ [0, T ], where the constant c σ is independent of s and t. Consequently, we have a continuous transition from the fractional Stefan-type problem to the classical Stefan-type problem as s 1 in the following sense.

Asymptotic Behaviour as t → ∞
In this section, we derive the asymptotic behaviour of the weak solutions as t → ∞, following the approach of the classical case in [20]. We first begin with a well-known asymptotic convergence result for the solutions of differential equations with maximal monotone operators.
Proposition 7 (See, for instance, Theorem 3.11 of [13]). Let ϕ be a lower semi-continuous convex functional on a Hilbert space H. Suppose that for all C ∈ R, the set {x ∈ H : . With this proposition, we can directly obtain the convergence of the generalised enthalpy solutions η(t) → η ∞ in the case whereg(t) =g ∞ for all t ≥ t 0 , i.e. the Dirichlet data is independent of time, with f (t) − f ∞ ∈ L 1 (t 0 , ∞; H −s (Ω)). For more generalg(t) converging to someg ∞ , we may also have a characterisation of the asymptotic behaviour of the generalised enthalpy solution towards the stationary solution, which can be written in terms of the stationary Dirichlet problem ϑ ∞ = g ∞ in Ω c for the temperature ϑ ∞ : Theorem 4. Let f ,g and η 0 satisfy the assumptions in Theorem 1 such that f − f ∞ ∈ L 1 (0, ∞; H −s (Ω)) ∩ L 2 (0, ∞; H −s (Ω)) andg −g ∞ ∈ W 1,1 (0, ∞; L 2 (Ω)) for given f ∞ ∈ H −s (Ω) andg ∞ ∈ H s (R d ). (We can subsequently define g ∞ and g(t) in the same spaces using (1.8) as explained in the Appendix A.) Let η ∈ β(ϑ) be the generalised enthalpy solution to the fractional Stefan-type problem (1.12) for all T > 0. Then, there exists an η ∞ ∈ L 2 (Ω) such that η(t) → η ∞ strongly in H −s (Ω) and weakly in L 2 (Ω) as t → ∞, Proof. We first note that, while η ∞ is not unique in general, there exists a unique weak temperature solution ϑ ∞ = γ(η ∞ ) to (5.1) with ϑ ∞ = g ∞ in Ω c by the Riesz representation theorem for A coercive and bounded, since we have the equivalent norms (1.5) in H s 0 (Ω). Furthermore, under our assumptions, by a similar approach to the Proposition 3.2 and its Corollary in [20], there is a positive constant M such that Let be any positive number. Since g is bounded, we can take a number t such that Also, let w be the solution of the fractional Stefan-type problem (1.12) corresponding Therefore, there is a number t ≥ t such that Also, as in Remark 5, we have that for some constant K for a.e. τ ≥ t . Integrating both sides over [t , t], we have for any t ≥ t . Therefore, if t, s ≥ t , This implies that η(t) converges in H −s (Ω) as t → ∞ to some η ∞ ∈ H −s (Ω). Also, since (5.5) holds for all t ≥ t and lim t→∞ w (t) = w ∞ , we have that w ∞ → η ∞ in H −s (Ω) as 0. Since w ∞ satisfies (5.4), so does η ∞ .
Next, taking a sequence {t n } n with t n → ∞ so that which is always possible by (5.6), we have, recalling that ϑ ∞ is the weak temperature solution to (5.1), Therefore, by (5.9), (5.13) and (5.12), Finally, since the duality in the left hand side of (5.12) is equivalent to the square of the H s 0 (Ω) norm of ϑ(t) − g(t) by (2.1), we may conclude the strong convergence result Remark 10. Similar asymptotic results as t → ∞ for the case s = 1 have been obtained in [20] considering other variants on the asymptotic behaviour of f andg. Earlier asymptotic behaviour results for s = 1 were obtained in Remarks 9 and 11 of [51] in the variational inequality form in a special case.

From Two Phases to One Phase
Let ν be a parameter such that (1.14) written with the Lipschitz graph γ ν corresponds to the two-phase problem when ν > 0, and to the one-phase problem when ν = 0. In this section, we obtain the solution to the one-phase problem, making use of the solution to the two-phase problem.
Consider the one-phase problem given with data Proof. We construct η o and ϑ o as the limit of an approximating sequence of η ν and ϑ ν . (See also the proof of Theorem A.1 in [24].) Indeed, since γ o is non-negative, Then, consider the strictly increasing approximation for ν > 0. Assuming γ o is Lipschitz continuous, so is γ ν . Also, γ ν clearly converges to γ o uniformly on compact sets as ν tends to zero. Furthermore, 3) is satisfied. The corresponding maximal monotone graph β ν is then given by β ν (r) = 1 ν (r − (Id + νβ o ) −1 (r)), (6.4) which is Lipschitz continuous with constant 1 ν . Therefore, from Theorem 1 and Theorem 2, we obtain the unique generalised enthalpy and weak temperature solutions η ν and ϑ ν of the approximate regularized problem with approximating compatible functions f ν , g ν and η ν 0 in the same spaces as the ones of the data such that η ν = β ν (ϑ ν ) are uniformly bounded in H 1 (0, T ; H −s (Ω) ∩ L ∞ (0, T ; L 2 (Ω)) for ν < 1, since the estimates (4.6)-(4.7) are independent of ν with for uniformly bounded η ν 0 , g ν (0) ∈ L 2 (Ω). We recall that (F n ) n∈N is an increasing set of finite dimensional subspaces of H s 0 (Ω), F * n = L(F n ) ⊂ H −s (Ω), and I F * n is the indicator function of F * n , i.e. I F * n = 0 in F * n , I F * n = +∞ elsewhere.

Remark 12.
Similarly to the convergence of the two-phase problem, it is possible to extend the results of Sections 4 and 5 to the one-phase problem.
with the Dirichlet boundary condition given by withg(t) defined on H s (R d ). Wheng ∈ BV (0, T ; H s (R d )) or H k (0, T ; H s (R d )) for k = 1, 2, by solving this Dirichlet problem, g will have the same time regularity asg.
Indeed, consider u = g −g. Then u satisfies u(t) = 0 in Ω c and Therefore, for any > 0, take a small enough h > 0 such that g(t)−g(t+h) < a * a * . Since is arbitrary, a.e. t in H s 0 (Ω), and the limit of the difference quotient is, by definition, ∂u ∂t . Therefore, ∂g ∂t = w(t) + ∂g ∂t (t), and we have that g has the same regularity asg in H 1 (0, T ; H s (R d )). Repeating this argument again by taking a second time derivative, we have the same result for g ifg ∈ H 2 (0, T ; H s (R d )).

B Appendix -The Variational Inequality Formulations
We observe that the formulation given in (1.14) can be formally transformed into a variational inequality formulation with fractional derivatives (see for example [41] or Chapter VII of [18]). Indeed, consider an element w ∈ H s 0 (Ω) independent of t and taking in (1.12) the test function ξ(x, τ ) = w(x) for τ ∈]t − , t + [ and ξ(x, τ ) = 0, dividing by 2 and letting → 0, denoting now by ·, · the duality between H −s (Ω) and H s 0 (Ω), we obtain Then, integrating with respect to time and using the regularity of η and its initial condition, we have, for almost all t ∈ [0, T ] and w ∈ H s 0 (Ω) by recalling that´t 0´R d AD s g · D s w = 0 for all w. We write η(t) = b(ϑ(t)) + λχ(t) for a.e. t for λ > 0 and b a given continuous and increasing function (see Figure 1). Then, denoting we observe that b(ϑ(t)) = b ∂Θ ∂t (t) ∈ L 2 (Ω) a.e. t. On the other hand, since H(r) is the subdifferential of the convex function r + , we have the inequality So, we obtain from (B.1) the nonlocal variational inequalitŷ for all w ∈ H s 0 (Ω) for a.e. t. By Theorem 1, ϑ − g ∈ L 2 (0, T ; H s 0 (Ω)), so Θ satisfies 3) with w =w(t) − ∂Θ ∂t (t), wherew(t) ∈ K(t), we obtain, for almost every t, which corresponds to the variational inequality formulations of Duvaut and Frémond (see [18], [50], [51] and [41]). With the same assumptions on f ,g and η 0 , we can obtain a solution Θ to (B.5), (B.4) using the Faedo-Galerkin method (refer to [51] or Chapter 3 of [41] for a proof starting from the variational inequality formulation (B.5), using the special basis of Appendix C. A similar result can also be obtained using the Rothe method (refer to Section 3.1 of [52]).
Similarly, for the one phase problem we can also obtain an equivalent variational inequality formulation, now of obstacle type. Indeed, governed by γ o , the weak temperature solution ϑ o obtained in (1.14 1ph ) is non-negative at all times t ∈ [0, T ]. Therefore, its primitive is also always non-negative, and satisfies when v(t), Θ o (t) ≥ 0. Therefore, we can rewrite the equation (B.1 o ) with w = v − Θ o (t) for v ∈ K + (t) as a variational inequality to obtain the following evolutionary obstacle-type problem for Θ o (t) ∈ K + (t): This corresponds to the nonlocal version of the parabolic variational inequality obtained by Duvaut [26] for the one-phase Stefan problem for the classical case s = 1. See also [40], [41] or [52].
C Appendix -Dependence of Eigenfunctions of L s A on 0 < s ≤ 1 Here we show the continuity of the eigenfunctions of L s A with respect to the parameter s, 0 < s ≤ 1. A similar result on s 1 can be found in Theorem 1.2 of [10] for the nonlocal p-Laplacian and Theorem 3. Then, by the Poincaré inequality, we have Therefore, for σ < s, By the estimate (C.2), for σ ≤ s → r ≤ 1, u s converges strongly to some u * in L 2 (Ω). As argued in Section 3.2 of [36], D s u s L 2 (R d ) d ≤ C for some constant C independent of s. Therefore, But by the a priori estimate on D s u s , ˆR This means that D r u * ∈ L 2 (R d ) d , and since Ω has a Lipschitz boundary, u * ∈ H r 0 (Ω). Furthermore, since D s · Φ → D r · Φ strongly in L 2 (R d ) d as s → r, sô Taking test functions ϕ ∈ C ∞ c (Ω), Extending this by density to all test functions v ∈ H r 0 (Ω), by the uniqueness of the solution to the homogeneous Dirichlet boundary problem (C.1) with s = r ≤ 1, we have that u * = u r . Therefore, for every h ∈ L 2 (Ω), T s (h) converges to T r (h) in L 2 (Ω) as s → r.
Theorem 8. Let 0 < σ ≤ s, r ≤ 1. For the sequence of operators T s : L 2 (Ω) → L 2 (Ω) given above, T s converges to T r strongly in the operator norm as s → r.
Proof. We first claim that, for each fixed s, it is possible to find an h s in the unit ball of L 2 (Ω) achieving the supremum, i.e. Indeed, for any maximizing sequence {h m } m , we can extract a subsequence which converges weakly to some h s which also belongs to the unit ball of L 2 (Ω). Since the embedding from L 2 (Ω) into H −σ (Ω) ⊂ H −s (Ω) ∩ H −r (Ω) is compact, and since T s and T r can also be considered continuous operators from H −s (Ω) into H s 0 (Ω) and from H −r (Ω) into H r 0 (Ω), respectively, both operators are also completely-continuous operators in L 2 (Ω), and so taking m to infinity we have the conclusion.
Having obtained the sequence {h s } s , since they are the weak limits of a uniformly bounded sequences, there exists h in the unit ball of L 2 (Ω) such that h s converge weakly in L 2 (Ω) and strongly in H −σ (Ω) to h. Then, by Lemma 3, for σ ≤ s, we have u H σ 0 (Ω) ≤ c σ u H s 0 (Ω) for u ∈ H s 0 (Ω) and consequently for the operator norm · s as an operator from H −s (Ω) to H s (Ω). Therefore, it follows that Similarly, we have Since T s (h) converges to T r (h) in L 2 (Ω) for every h ∈ L 2 (Ω), for any > 0, we can pick a δ > 0 such that, for |s − r| ≤ δ, we have h s − h H −σ (Ω) ≤ a * 4c 2 σ and T s (h) − T r (h) L 2 (Ω) ≤ 2 . Therefore, As a corollary, by Theorem 2.3.1 of [29], we have Corollary 2. For the operators T s , T r as given in the previous theorem, let λ s k = λ s k (T s ) and λ r k = λ r k (T r )be the k-th eigenvalues of T s and of T r respectively for s and for r, 0 < σ ≤ s, r ≤ 1. Then, |λ s k − λ r k | ≤ T s − T r := sup f L 2 (Ω) ≤1 (T s − T r )(f ) .