Subordinated Gaussian Random Fields in Elliptic Partial Differential Equations

To model subsurface flow in uncertain heterogeneous\ fractured media an elliptic equation with a discontinuous stochastic diffusion coefficient - also called random field - may be used. In case of a one-dimensional parameter space, L\'evy processes allow for jumps and display great flexibility in the distributions used. However, in various situations (e.g. microstructure modeling), a one-dimensional parameter space is not sufficient. Classical extensions of L\'evy processes on two parameter dimensions suffer from the fact that they do not allow for spatial discontinuities. In this paper a new subordination approach is employed to generate L\'evy-type discontinuous random fields on a two-dimensional spatial parameter domain. Existence and uniqueness of a (pathwise) solution to a general elliptic partial differential equation is proved and an approximation theory for the diffusion coefficient and the corresponding solution provided. Further, numerical examples using a Monte Carlo approach on a Finite Element discretization validate our theoretical results.

1. Introduction. Over the last decade partial differential equations with stochastic opera-tors\ data\ domain became a widely studied object. This branch of research is oftentimes called uncertainty quantification. Especially for problems where data is sparse or measurement errors are unavoidable, like subsurface flow problems, the theory provides an approach to quantify this uncertainty. There are two main approaches to discretize the uncertain problem: intrusive and non-intrusive methods. The former require the solution of a high dimensional partial differential equation, where the dimensionality depends on the smoothness of the random field or process (see for example [8], [21], [27] and the references therein). The latter consist of (essentially) sampling methods and require multiple solutions of a low dimensional problem (see, among others, [1], [10], [12], [13], [26], [31]). Up to date mainly Gaussian random fields were used to model the diffusivity in an elliptic equation (as a model for a subsurface flow problem). Gaussian random fields have the advantage that they may be used in both approaches and that they are stochastically very well understood objects. A great disadvantage is however, that the distributions underlying the field are Gaussian and therefore the field lack flexibility, in the sense that the field is continuous and cannot have pointwise marginal distributions with heavy-tails.
In this paper we propose a two-dimensional subordinated Gaussian random field as stochastic diffusion coefficient in an elliptic equation. The subordinated Gaussian random field is a type of a (discontinuous) Lévy field. Different subordinators display unique patterns in the discontinuities and have varied marginal distributions (see [9]). Naturally the spatial regularity of a subordinated Gaussian random field depends on the subordinator. We prove existence and uniqueness of a solution to the elliptic equation in a pathwise sense and provide different discretization schemes.
We structured the rest of the paper as follows: In Section 2 we introduce a general pathwise existence and uniqueness result for a stochastic elliptic equation under mild assumptions on the coefficient. These assumptions accommodate the subordinated Gaussian random fields we introduce in Section 3. In Section 4 we approximate the specific diffusion coefficient which is used in this paper and show convergence of the elliptic equation with the approximated coefficient to the unapproximated solution in Section 5. Section 6 provides spatial approximation methods and in Section 7 numerical examples are presented.
2. The stochastic elliptic problem. In this section we introduce the framework of the general stochastic elliptic boundary value problem which allows for discontinuous diffusion coefficients. For the general setting and pathwise existence theory we follow [12]. In the following, let (Ω, F, P) be a complete probability space. To accommodate Banach-space-valued random variables we introduce the so-called Bochner spaces. Definition 2.1. Let (B, ⋅ B ) be a Banach space and Z be a B valued random variable, i.e. a strongly measurable function Z ∶ Ω → B. The space L p (Ω; B) contains all B-valued random variables with Z L p (Ω;B) < +∞, for p ∈ [1, +∞), where the norm is defined by 2.1. Problem formulation. Let D ⊂ R d for d ∈ N be a bounded, connected Lipschitz domain. We consider the equation −∇(a(ω, x)∇u(ω, x)) = f (ω, x) in Ω × D, (2.1) where a ∶ Ω×D → R + is a stochastic (jump diffusion) coefficient and f ∶ Ω×D → R is a (measurable) random source function. Further, we impose the following boundary conditions where we assume to have a decomposition ∂D = Γ 1 . ∪ Γ 2 with two (d − 1)-dimensional manifolds Γ 1 , Γ 2 such that the exterior normal derivative → n ⋅ ∇u on Γ 2 is well-defined for every u ∈ C 1 (D).
Here, → n is the outward unit normal vector to Γ 2 and g ∶ Ω × Γ 2 → R a measurable function. Note that we just reduce the theoretical analysis to the case of homogeneous Dirichlet boundary conditions to simplify notation. It would be also possible to work under non-homogeneous Dirichlet boundary conditions, since such a problem can always be considered as a version of (2.1) - (2.3) where the source term and the Neumann data have been changed (see also [12,Remark 2.1]).
We now state assumptions under which the elliptic boundary value problem has a unique solution.
Remark 2.3. Note that Assumption 2.2 implies that the real-valued mappings a − , a + ∶ Ω → R are measurable. This can be seen as follows: For fixed p ≥ 1 consider the mapping which is well-defined by Assumption 2.2. It follows from the definition of the Lebesgue integral and Assumption 2.2 i that the mapping ω ↦ I p (ω) is F − B(R) measurable. For a fixed ω ∈ Ω, by the embedding theorem for L p spaces (see [2,Theorem 2.14]), we get a + (ω) = lim m→∞ a(ω, ⋅) L m (D) .
Since this holds for all ω ∈ Ω we obtain by [4,Lemma 4.29] that the mapping ω ↦ a + (ω), is F − B(R)-measurable. The measurability of ω ↦ a − (ω) follows analogously. Note that we do not treat the random coefficient a ∶ Ω×D → R as a L ∞ (D)-valued random variable, since L ∞ (D) is not separable and therefore the strong measurability of the mapping a ∶ Ω → L ∞ (D) is only guaranteed in a very restrictive setting. Nevertheless, the measurability of the functions a + and a − allows taking expectations of these real-valued random variables. In order to avoid confusion about that, we use the notation E(ess sup x∈D ⋅ s ) 1 s .

Weak solution.
We denote by H 1 (D) the Sobolev space on D with the norm (see for example [20,Section 5.2]). Here, 2 denotes the Euclidean norm of the vector x ∈ R d . Further, we denote by T the trace operator with [18]). We define the subspace V ⊂ H 1 (D) as with the standard Sobolev norm, i.e. ⋅ V ∶= ⋅ H 1 (D) .
We multiply Equation (2.1) by a test function v ∈ V , integrate by parts and use the boundary conditions (2.2) and (2.3) to obtain This leads to the following pathwise weak formulation of the problem: For any ω ∈ Ω, given f (ω, ⋅) ∈ V ′ and g(ω, for all v ∈ V . The function u(ω, ⋅) is then called pathwise weak solution to problem (2.1) -(2.3). Here, the bilinear form B a(ω) and the operator F ω are given by for fixed ω ∈ Ω, where the integrals in F ω are understood as the duality pairings: Theorem 2.4. (see [12,Theorem 2.5]) Under Assumption 2.2, there exists a unique pathwise weak solution u(ω, ⋅) ∈ V to problem (2.4) for very ω ∈ Ω. Furthermore, u ∈ L r (Ω; V ) and u L r (Ω;V ) ≤ C(a − , D, p)( f L q (Ω;H) + g L q (Ω;L 2 (Γ2)) ), where C(a − , D, p) > 0 is a constant depending on a − , p and the volume of D.
In addition to the existence of the solution, the following remark gives a rigorous justification for the measurability of the solution mapping which maps any ω ∈ Ω on the corresponding pathwise weak PDE solution.
Remark 2.5. Let (v n , n ∈ N) ⊂ V be an orthonormal basis of the separable Hilbert space V . For every n ∈ N we define the mapping It is easy to see that this mapping is Carathéodory for any n ∈ N, i.e. J n (⋅, v) is F-B(R)-measurable for any fixed v ∈ V and J n (ω, ⋅) is continuous on V for any fixed ω ∈ Ω. We define the correspondences for every n ∈ N. It follows from [4,Corollary 18.8] that this correspondence has a measurable graph, i.e.
Further, by Assumption 2.2 and the Lax-Milgram Lemma (see for example [25,Lemma 6.97] and [12, Theorem 2.5]) we know that for every fixed ω ∈ Ω, there exists a unique solution u(ω, ⋅) ∈ V satisfying (2.4) for every v ∈ V . Therefore, we obtain for the graph of the solution mapping: This implies for an arbitrary measurable setṼ ∈ B(V ) and therefore {ω ∈ Ω u(ω, ⋅) ∈Ṽ } ∈ F, by the projection theorem (see [4,Theorem 18.25]), which gives the measurability of the solution mapping.
3. Subordinated Gaussian random fields. A random field W ∶ Ω × D → R is called R − valued Gaussian random field (GRF) if for any tuple (x 1 , . . . , x n ) ⊂ D and any number n ∈ N the R n -valued random variable is multivariate normally distributed (see [3,Section 1.2]). Here x T denotes the transpose of the vector x. We denote by the associated mean and covariance function. The covariance operator Q ∶ L 2 (D) → L 2 (D) of W is defined by Further, if D ⊂ R d is compact and W is centered, i.e. m ≡ 0, there exists a decreasing sequence (λ i , i ∈ N) of real eigenvalues of Q with corresponding eigenfunctions (e i , i ∈ N) ⊂ L 2 (D) which form an orthonormal basis of L 2 (D) (see [3,Section 3.2] and [32, Theorem VI.3.2 and Chapter II.3]).
Remark 3.2. Let B = (B(t), t ≥ 0) be a standard Brownian motion and l = (l(t), t ≥ 0) be a Lévy subordinator. The stochastic process defined by is called subordinated Brownian motion and is again a Lévy process (see [5,Theorem 1.3.25]).
In [9] the authors propose a new approach to extend standard subordinated Lévy processes on a higher dimensional parameter space. Motivated by the rich class of subordinated Brownian motions, the authors construct discontinuous random fields by subordinating a GRF on a d-dimensional parameter domain by d one-dimensional Lévy subordinators. In case of a twodimensional parameter space the construction is as follows: For a GRF W ∶ Ω × R 2 + → R and two (Lévy-)subordinators l 1 , l 2 on [0, D], with a finite D > 0, we define the real-valued random field L(x, y) ∶= W (l 1 (x), l 2 (y)), for x, y ∈ [0, D]. This construction yields a rich class of discontinuous random fields which also admit a Lévy-Khinchin-type formula. Further, the newly constructed random fields are also interesting for practical reasons, since they admit a semi-explicit formula for the covariance function which is very useful for applications, e.g. in statistical fitting. For a theoretical investigation of the constructed random fields we refer to [9]. 3.2. Subordinated GRFs as diffusion coefficients in elliptic problems. In the following, we define the specific diffusion coefficient that we consider in problem (2.1) -(2.3). In order to allow discontinuities, we incorporate a subordinated GRF in the coefficient additionally to a Gaussian component. The construction of the coefficient is done so that Theorem 2.4 is applicable and, at the same time, the coefficient is as versatile as possible. which are independent of the GRFs W 1 and W 2 .
Remark 3.4. The first two assumptions ensure that the diffusion coefficient a is positive over the domain D. To show the convergence of the approximated diffusion coefficient in Subsection 5.1 we have to impose independence of the GRFs W 1 and W 2 (see Assumption 5.1 and the proof of Theorem 5.4). This assumption is in the sense natural as also one-dimensional Lévy processes admit an additive decomposition into a continuous part and a pure-jump part which are stochastically independent (Lévy-Itô decomposition, see e.g. [5,Theorem 2.4.11]). For the same reason the assumption that the Lévy subordinators are independent of the GRFs is also natural (see for example [5,Section 1.3.2]).
In order to verify Assumption 2.2 i and ii we need the following Lemma.
Theorem 3.6. Let a be as in Definition 3.3 and let f ∈ L q (Ω; H), g ∈ L q (Ω; L 2 (Γ 2 )) for some q ∈ [1, +∞). Then there exists a unique pathwise weak solution u(ω, ⋅) ∈ V to problem (2.1) for P-almost every ω ∈ Ω. Furthermore, u ∈ L r (Ω; V ) for all r ∈ [1, q) and where C(a − , D) > 0 is a constant depending only on the indicated parameter and the volume of D.

Approximation of the diffusion coefficient.
To simulate the solution to the elliptic equation we need to define a tractable approximation of the diffusion coefficient.
In order to approximate the solution to problem (2.1) -(2.3) we face a new challenge regarding the GRF W 2 which is subordinated by the Lévy processes l 1 and l 2 : Due to the fact that the Lévy subordinators in general can attain any value in [0, +∞) we have to consider (and approximate) the GRF W 2 on the unbounded domain [0, +∞). In most cases where elliptic PDEs of the form (2.1) have been considered with a random coefficient, the problem is stated on a bounded domain, see e. g. [12,15,14,22,24]. Many regularity results for GRFs formulated for a bounded parameter space cannot easily be transferred to an unbounded parameter space (see also [3,Chapter 1], especially the discussion on p. 13). Even the Karhunen-Loève expansion of a GRF requires compactness of the domain (see e.g. [3, Section 3.2]).
Furthermore, to show convergence of the solution in Section 5 we need to bound the coefficient from above by a deterministic upper bound A (see Theorem 5.9 Remark 5.7). Subsequently we show that this induces an error in the solution approximation which can be controlled and which vanishes for growing A (see Subsection 5.1).
Remark 4.1. We note that the influence of this problem modification can be controlled: one may choose K > 0 such that for any ε > 0. In other words, pathwise the modified problem coincides with the original one up to a set of samples, whose probability can be made arbitrarily small.

4.2.
Second modification: diffusion cut-off. We consider again the cut function χ A (z) ∶= min(z, A) for z ∈ [0, +∞) with a fixed positive number A > 0 and consider the following problem where we impose the boundary conditions The diffusion coefficient a K,A is defined by Again, Theorem 3.6 applies in this case and yields the existence of a pathwise weak solution u K,A ∈ L r (Ω; V ) for r ∈ [1, q) if f ∈ L q (Ω; H) and g ∈ L q (Ω; L 2 (Γ 2 )). The error of the modification vanishes for growing A as shown in the end of Section 5.1.

Approximation of GRF and subordinators.
Here, we show how to approximate the modificated diffusion coefficient a K,A using approximations W ε W 1 ≈ W 1 , W ε W 2 ≈ W 2 of the GRFs and l ε l 1 ≈ l 1 , l ε l 2 ≈ l 2 of the Lévy subordinators. To this end we need additional assumptions on the data of the elliptic problem, the covariance operators of the Gaussian fields and the subordinators.
the covariance functions of these random fields and by Q 1 , Q 2 the associated covariance operators defined by i , i ∈ N) the eigenpairs associated to the covariance operators Q 1 and Q 2 . In particular, (e i We assume that the eigenfunctions are continuously differentiable and there exist positive constants α, β, C e , C λ > 0 such that for any i ∈ N it holds ii There exist constants φ, ψ, C lip > 0 such that the continuous functions iii f ∈ L q (Ω; H) and g ∈ L q (Ω; L 2 (Γ 2 )) for some q ∈ (1, +∞). iv a ∶ D → (0, +∞) is deterministic, continuous and there exist constants a + , a − > 0 with a − ≤ a(x, y) ≤ a + for (x, y) ∈ D. v l 1 and l 2 are Lévy subordinatos on [0, D] with Lévy triplets (γ 1 , 0, ν 1 ) and (γ 2 , 0, ν 2 ) which are intependent of the GRFs W 1 and W 2 . Further, we assume that we have approximations l of these processes and there exist constants η, C l > 0 such that for every s ∈ [1, η − 1) it holds Remark 4.3. Note that the first assumption on the eigenpairs of the GRFs is natural (see [12] and [22]). For example, the case that Q 1 , Q 2 are Matérn covariance operators are included. Assumption 4.2 ii is necessary to be able to quantify the error of the approximation of the diffusion coefficient. Assumption 4.2 iii is necessary to ensure the existence of a solution and has already been formulated in Assumption 2.2. The last assumption ensures that we can approximate the Lévy subordinators in an L s -sense. This can always be achieved under appropriate assumptions on the tails of the distribution of the subordinators, see [11, Assumption 3.6, Assumption 3.7 and Theorem 3.21].
For any numerical simulation we have to approximate the GRF as well as the subordinating Lévy processes, which results in an additional approximation of the coefficient a K,A given in Equation (4.8). In the following we want to quantify the error induced by this approximation.
Next, we prove a bound on the error of the approximated diffusion coefficient, where the GRFs are approximated by a discrete evaluation and (bi-)linear interpolation between these points (see [23] and [24]).
which are constructed by point evaluation of the random fields W 1 and W 2 on the grids and linear interpolation between the grid points. Under Assumption 4.2 i it holds for n ∈ [1, +∞): where β and α are the parameters from Assumption 4.2.
Proof. Note that for any fixed ω ∈ Ω and is constructed by (bi-)linear interpolation of the GRF W 1 and the piecewise linear interpolants attain their maximum and minimum at the corners (the Hessian evaluated at the (unique) stationary point of the bilinear basis functions is always indefinite). Therefore, for a fixed (x, y) Using this observation we estimate where we used Equation (4.9) in the last step. Equivalently the error bound for W 2 follows.

Convergence to the modificated diffusion coefficient. Given some approximations W
≈ l 2 as in Assumption 4.2 v as well as some fixed constants K, A > 0, we approximate the diffusion coefficient for (x, y) ∈ D. To prove a convergence result for this approximated coefficient (Theorem 4.8) we need the following two technical lemmas. The second can be proved by the use of [29,Proposition 1.16]. For a detailed proof we refer to [9].
Proof. The case n = m is trivial. For n > m we use Hölder's inequality and obtain s. continuous random field and let Z ∶ Ω → R d + be a R d + -valued random variable which is independent of the random field W . Further, let ϕ ∶ R → R be a deterministic, continuous function. It holds ≈ W 2 be approximations of the GRFs on discrete grids as in Lemma 4.4. Further, let 1 ≤ t ≤ s < η − 1 and 0 < γ < min(1, β (2α)) such that sγ ≥ 2. Under Assumption 4.2 we get the following error bound for the approximation of the diffusion coefficient: with a constant C which does not depend on the discretization parameters ε W and ε l .
Proof. Since the cut function χ A is Lipschitz continuous with Lipschitz constant 1 we calculate 2 ))) L s (Ω;L t ([0,D] 2 )) =∶ I 1 + I 2 First, we consider I 1 and use Assumption 4.2 ii and the same calculation as in Remark 4.5 to get The mean value theorem yields, for fixed (x, y) ∈ [0, D] 2 and an appropriately chosen value ξ ∈ (min(W 1 (x, y), W for P-almost every ω ∈ Ω. As already mentioned in the proof of Lemma 4.4, for any ω ∈ Ω and Therefore, we obtain the pathwise estimate Finally, we obtain for any n 1 , n 2 ∈ [1, +∞) with 1 n 1 + 1 n 2 = 1 by Hölder's inequality where we used Lemma 4.4 and the fact that We use the same calculation as in Remark 4.5 and the Lipschitz continuity of Φ 2 to calculate for the summand I 4 where we used Lemma 4.4. It remains to bound the summand I 3 : We estimate using Lemma 4.6 by Equation (4.9). Further, we know from Hölder's inequality for γs ≥ 2 that it holds and therefore we calculate where we used the Lipschitz continuity of χ K and Assumption 4.2 v in the last step. Therefore, we finally obtain 5. Convergence analysis. In this section we derive an error bound for the approximation of the solution. We split the error in two components: the first component is associated with the cut-off of the diffusion coefficient we described in Subsection 4.2. The second error contributor corresponds to the approximation of the GRFs and the Lévy subordinators we considered in Subsection 4.4.
Let r ∈ [1, q) with q as in Assumption 4.2 iii and denote by u K ∈ L r (Ω; V ) the weak solution to problem (4.1) -(4.4). Further, let u Note that Theorem 3.6 also applies to the elliptic problem with coefficient a . The aim of this section is to quantify the error of the approximation u By the triangle inequality we obtain for an arbitrary norm ⋅ (to be specified later). Here, u K,A is the solution to the truncated problem (4.5) -(4.8). We consider the two error contributions E 1 and E 2 separately.

Bound on E 1 .
Assumption 5.1. We assume that the GRFs W 1 and W 2 occurring in the diffusion coefficient (3.1) are stochastically independent.
The aim of this subsection is to show that the first error contributor E 1 in Equation (5.4) vanishes for increasing cut-off threshold A. The strategy consists of two separated steps: in the first step we show the stability of the solution, which means that the value E 1 can be controlled by the quality of the approximation of the diffusion coefficient a K,A ≈ a K . In the second step, we show that the quality of the approximation of the diffusion coefficient can be controlled by the cut-off threshold A. The first step is given by Theorem 5.3. In order to prove it we need the following lemma.
Lemma 5.2. For fixed cut-off levels A and K we consider the solution u K ∈ L r (Ω; V ) and its approximation u K,A ∈ L r (Ω; V ) for r ∈ [1, q). It holds the pathwise estimate for P-almost every ω ∈ Ω. Here, the constant C(a − , D) only depends on the indicated parameters and we define a K,+ (ω) ∶= max{1, ess sup (x,y)∈D a K (ω, x, y)} < ∞ for ω ∈ Ω.
Proof. For a fixed ω ∈ Ω we consider the variational problem: find a uniqueŵ ∈ V such that for all v ∈ V . By the Lax-Milgram theorem there exists a unique solutionŵ ∈ V with (see Theorem 3.6 and [12, Theorem 2.5]). Therefore, we obtain by Hölder's inequality where we used Young's inequality in the last step. Finally, we obtain Theorem 5.3. Let f ∈ L q (Ω; H) and g ∈ L q (Ω; L 2 (Γ 2 )) for some q ∈ [1, +∞). Further, for a given number t ∈ (1, +∞) we define the dual number n ∶= t t−1 . Then, for any for r ∈ [1, q n), holds Proof. By a direct calculation we obtain We estimate using Hölder's inequality and therefore we obtain Using Lemma 5.2 we obtain the pathwise estimate Using again Hölder's inequality we have By assumption it holds nr < q. Therefore, we can choose a real number ρ > 1 such that nrρ < q. We define the dual number ρ ′ ∶= ρ ρ−1 ∈ (1, +∞) and use Hölder's inequality to calculate Obviously Theorem 3.6 applies also to problem (4.1) -(4.4). Therefore, since nrρ < q by assumption we conclude where we additionally used the fact that E(a nrρ ′ K,+ ) 1 nrρ ′ < +∞ (see Theorem 5.4).
In other words, finding a bound for the error contribution E 1 in Equation (5.4) reduces to quantifying the quality of the approximation of the diffusion coefficient a K,A ≈ a K . For readability the proof of the following theorem can be found in Appendix A Theorem 5.4. For any n ∈ (1, +∞) it holds Further, for any δ > 0 there exists a constant A = A(δ, n) > 0 such that

Bound on E 2 . The aim is to bound the second term of Equation (5.4) given by
in an appropriate norm. For technical reasons we have to impose an additional assumption on the solution of the truncated problem. The subsequent remarks discuss situations under which this assumption is fulfilled.
Assumption 5.5. We assume that there exist constants j reg > 0 and k reg ≥ 2 such that Since u K,A ∈ H 1 (D) we already know that ∇u K,A ∈ L 2 (D). Assumption 5.5 requires a slightly higher integrability over the spatial domain. Since this is an assumption on the regularity of the solution u K,A we denote the above constant by C reg .
Remark 5.7. Note that there are several results about higher integrability of the gradient of the solution to an elliptic PDE of the form (4.5) -(4.8). For instance [28] yields that the solution u K,A has H 1+δ (2π) regularity with δ = min(1, a − , A −1 ) under mixed boundary conditions and under the assumption that a K,A is piecewise constant (see [28,Theorem 7.3]). This corresponds to the case where no Gaussian noise is considered (i.e. Φ 1 ≡ 0) and a is constant. Another important result is given in [19]. It follows by [19,Theorem 1] that under the assumption that there exists q > 2 with f ∈ L q (D) P − a.s. there exists a constant C = C(D, f L q (D) , a − , A) and a positive number ϑ = ϑ(D, f L q (D) , a − , A) > 0 only depending on the indicated parameters, such that: In particular, if the right hand side f of the problem is deterministic, then ϑ and the constant C in (5.7) are deterministic and one immediately obtains for any k reg ≥ 1 and a deterministic, positive constant ϑ > 0. We note that although [19,Theorem 1] suggests a dependence of the parameter ϑ = ϑ(D, f L q (D) , a − , A) on the constant A, this dependence is not numerically detectable for our diffusion coefficient, as numerical experiments show. Of course, it depends on the other parameters D, f L q (D) and a − .
Next, we show that for a given approximation of the diffusion coefficient, the resulting error contributor E 2 is bounded by the approximation error of the diffusion coefficient. Similar to the corresponding assertion we gave in Subsection 5.1 we need the following Lemma for the proof of this error bound. For a proof we refer to Lemma 5.2.
Theorem 5.9. Let r ≥ 2 and b, c ∈ [1, +∞] be given such that it holds rcγ ≥ 2 and 2b ≤ rc < η − 1 with a fixed real number γ ∈ (0, min(1, β (2α)). Here, the parameters η, α and β are determined by the GRFs W 1 , W 2 and the Lévy subordinators l 1 , l 2 (see Assumption 4.2). Let m, n ∈ [1, +∞] be real numbers such that and let k reg ≥ 2 and j reg > 0 be the regularity specifiers given by Assumption 5.5. If it holds that n < 1 + j reg 2 and rm < k reg , then the approximated solution u (ε W ,ε l ) K,A converges to the solution u K,A of the truncated problem for ε W , ε l → 0 and it holds Proof. By a direct calculation we obtain the pathwise estimate ) is the weak solution to problem (4.5) -(4.8) (resp. (5.1) - Using Hölder's inequality we calculate and therefore Next, we apply Lemma 5.8 to obtain the following estimate.
Hence, it remains to bound the norm (a (ε W ,ε l ) K,A − a K,A )∇u K,A L 2 (D) . By Hölder's inequality we obtain Applying Hölder's inequality once more we estimate By assumption it holds rm < k reg . Hence, we can choose a real number ρ > 1 such that rmρ < k reg . We define the dual number ρ ′ ∶= ρ 1−ρ ∈ (1, +∞) and use Hölder's inequality to obtain where we again used the fact that E(a nrρ ′ K,+ ) 1 nrρ ′ < +∞ (see Theorem 5.4) together with Assumption 5.5. Finally we obtain the estimate , where we applied Theorem 4.8 in the last estimate.
We close this section with a remark on how to choose the parameters A, ε W and ε l to obtain an approximation error smaller than any given threshold δ > 0.
Remark 5.10. For any given parameter K large enough (see Remark 4.1), we choose a positive numer A > 0 such that the first error contributor satisfies E 1 = u K − u K,A L r (Ω;V ) < δ 2 (see Theorem 5.3 and Theorem 5.4). Afterwards, under the assumptions of Theorem 5.9, we may choose the approximation parameters ε W and ε l small enough, such that the secontd error contributor satisfies E 2 = u K,A − u (ε W ,ε l ) K,A L r (Ω;V ) < δ 2. Hence, we get an overall error smaller than δ (see Equation (5.4)). where the approximated diffusion coefficient a (ε W ,ε l ) K,A is given by (4.10). Therefore, for almost all ω ∈ Ω, we have to find a function u for every v ∈ V . Here, K, A, ε W , ε l are fixed approximation parameters. In order to solve this variational problem numerically we consider a standard Galerkin approach and assume V = (V , ∈ N 0 ) to be a sequence of finite-dimensional subspaces V ⊂ V with dim(V ) = d and V ⊂ V +1 for all ≥ 0. We denote by (h , ∈ N 0 ) the corresponding sequence of refinement sizes which is assumed to decrease monotoncally to zero for → ∞. Let ∈ N 0 be fixed and denote by {v The (pathwise) discrete version of (6.1) reads: We expand the function u where the coefficient vector c = (c 1 , . . . , c d ) T ∈ R d is determined by the linear equation system with a stochastic stiffness matrix B(ω) i,j = B a (ε W ,ε l ) Remark 6.1. Let (K , ∈ N 0 ) be a sequence of triangulations on D and denote by θ > 0 the minimum interior angle of all triangles in K . We assume θ ≥ θ > 0 for a positive constant θ and define the maximum diameter of the triangulation K by h ∶= max K∈K diam(K), for ∈ N 0 . Further, we define the finite dimensional subspaces by V ∶= {v ∈ V v K ∈ P 1 , K ∈ K }, where P 1 denotes the space of all polynomials up to degree one. If we assume that for P−almost all ω ∈ Ω it holds u (ε W ,ε l ) K,A (ω, ⋅) ∈ H 1+κa (D) for some positive number κ a > 0, the pathwise discretization error is bounded by Céa's lemma P-a.s. by (see [12,Section 4] and [25,Chapter 8]). If the bound u (ε W ,ε l ) K,A L 2 (Ω;H 1+κa (D)) ≤ C u = C u (K, A) is finite for the fixed approximation parameters K, A, we immediately obtain Note that, for general jump-diffusion problems, one obtains a discretization error of order κ a ∈ (1 2, 1). In general, we cannot expect the full order of convergence κ a = 1 since the diffusion coefficient is discontinuous. Without special treatment of the interfaces with respect to the triangulation, one cannot expect a convergence rate which is higher than κ a = 1 2 for the deterministic problem (see [7] and [12]).

Sample-adapted triangulations.
In [12], the authors suggest sample-adapted triangulations to improve the convergence rate of the FE approximation: Consider a fixed ω ∈ Ω and assume that the discontinuities of the diffusion coefficient are described by the partition T (ω) = (T i , i = 1, . . . , τ (ω)) of the domain D where τ (ω) describes the number of elements in the partition. We consider finite-dimensional subspacesV (ω) ⊂ V with (stochastic) dimension d (ω) ∈ N. We denote byθ (ω) the minimal interior angle within K (ω) and assume the existence of a positive number θ > 0 such that inf{θ (ω) ∈ N 0 } ≥ θ for P-almost all ω ∈ Ω. Assume that K (ω) is a triangulation of D which is adjusted to the partition T (ω) in the sense that for every i = 1, . . . , τ (ω) it holds is a deterministic, decreasing sequence of refinement thresholds which converges to zero (see Figure 2). Sample-adapted triangulations lead to an improved convergence rate for our problem (see [12,Section 4.1] and Section 7). This observation together with Remark 6.1 motivate the following assumption for the rest of this paper (see [12,Assumption 4.4]).
Assumption 6.2. There exist deterministic constantsĈ u , C u ,κ a , κ a > 0 such that for any ε W , ε l > 0 and any ∈ N 0 , the Finite Element approximation errors ofû where the constantsĈ u , C u may depend on a, f, g, K, A but are independent ofĥ , h ,κ a and κ a .
In practice, as we can see in the numerical examples in the subsequent section, one can (at least) recover the deterministic rates in the strong error. In fact, in the non-adapted case it is possible to get better convergence rates than expected for some examples. By construction of our random field, we always obtain an interface geometry with fixed angles and bounded jump height, which have great influence on the solution regularity, see e.g. [28].

Numerical examples.
In this section we verify our theoretical results in numerical examples. In all experiments we work on the domain D = (0, 1) 2 and use a FE method with hat-function basis. Here, we distinguish between the standard FEM approach and the sample-adapted FEM approach introduced in Section 6. We compare both and investigate how different Lévy subordinators influence the strong convergence rate. In our first example we use Poisson processes with low intensity to investigate the superiority of the presented sample-adapted triangulation. In the second example we use Poisson subordinators with a significantly higher intensity. Besides Poisson subordinators we also use Gamma processes which have infinite activity.
7.1. Strong error approximation. In each numerical experiment we choose a problem dependent cut-off level K for the subordinators in (4.4) large enough so that its influence is negligibly (see Remark 4.1). Further, we choose the cut-off level for the diffusion coefficient A in (4.8) large enough such that it has no influence in numerical experiments, e.g. A = 50 and therefore the error induced by the error contributor E 1 in (5.4) can be neglected in our experiments. We estimate the strong error using a standard Monte Carlo estimator. Assume that a sequence of (sample-adapted) finite-dimensional subspaces (V , ∈ N 0 ) ⊂ V is given where we use the notation of Section 6. For readability we only treat the case of pathwise sample-adapted Finite Element approximations in the rest of the theoretical consideration in this subsection. We would like to point out, however, that similar arguments lead to the corresponding results for standard FE approximations.
Under the assumptions of Theorem 5.9 and Assumption 6.2 we obtain with a constant C = C(C reg , D, a − ,Ĉ u ). Therefore, in order to equilibrate all error contributions, we choose the approximation parameters ε W and ε l in the following way: For readability, we omit the cut-off parameters K and A in the following and use the notation u ,ε W ,ε l =û (ε W ,ε l ) K,A, . Choosing the approximation parameters ε W , ε l according to (7.2), we can investigate the strong error convergence rate by a Monte Carlo estimation of the left hand side of (7.1): for a fixed natural number M ∈ N we approximate where (u for ω ∈ Ω. We choose W 1 to be a Matérn-1.5-GRF on D with correlation length r 1 = 0.5 and different variance parameters σ 2 1 . Further, we set W 2 to be a Matérn-1.5-GRF on [0, K] 2 which is independent of W 1 with different variances σ 2 2 and correlation lengths r 2 . We use a reference grid with 800 × 800 equally spaced points on the domain D for interpolation and prolongation.

Poisson subordinators.
In this section we use Poisson processes to subordinate the GRF W 2 in the diffusion coefficient in (3.1). We consider both, high and low intensity Poisson processes and vary the boundary conditions. Further, using Poisson subordinators allows for a detailed investigation of the approximation error caused by approximating the Lévy subordinators l 1 and l 2 according to Assumption 4.2 v. 7.3.1. The two approximation methods. Using Poisson processes as subordinators allows for two different simulation approaches in the numerical examples: the first approach is an exact and grid-independent simulation of a Poisson process using for example the Method of Exponential Spacings or the Uniform Method (see [30,Section 8.1.2]). On the other hand, one may also work with approximations of the Poisson processes satisfying Assumption 4.2 v.
We sample the values of the Poisson(λ)-processes l 1 and l 2 on an equidistant grid {x i , i = 0, ..., N l } with x 0 = 0 and x N l = 1 and step size x i+1 − x i ≤ ε l ≤ 1 for all i = 0, . . . , N l − 1. Further, we approximate the stochastic processes by a piecewise constant extension l (ε l ) j ≈ l j of the values on the grid: for j = 1, 2. Since the Poisson process has independent, Poisson distributed increments, values of the Poisson process at the discrete points {x i , i = 0, . . . , N l } can be generated by adding independent Poisson distributed random variables. In the following we refer to this approach as the approximation approach to simulate a Poisson process. Note that in this case Assumption 4.2 v holds with η = +∞. In fact, for any s ∈ [1, +∞) we obtain for j = 1, 2 and an arbitrary x which is independent of the specific x ∈ [0, 1). Note that this also holds for x = D = 1 and therefore For a Poisson process with parameter λ we obtain where the series converges by the ratio test.
Since the Poisson process allows for both approaches -approximation and exact simulation of the process -the use of these processes are suitable to investigate the additional error in the approximation of the PDE solution resulting from an approximation of the subordinators. 7.3.2. Poisson subordinators: low intensity and mixed boundary conditions. In this example we choose l 1 and l 2 to be Poisson(1)-subordinators. Further, the variance parameter of the GRF W 1 is set to be σ 1 = 1.5 and the variance and correlation parameters of the GRF W 2 are given by σ 2 = 0.3 and r 2 = 1.
We approximate the GRFs W 1 and W 2 by the circulant embedding method (see [23] and [24]) to obtain approximations W ≈ W 2 as in Lemma 4.4. Since η = +∞ and f ∈ L q (Ω; H) for every q ≥ 1 we choose for any positive δ > 0 r = 2, c = b = 1 + δ to obtain from Theorem 5.9 where we have to assume that j reg ≥ 2((1 + δ) δ − 1) and k reg ≥ 2(1 + δ) δ for the regularity constants j reg , k reg given in Assumption 5.5. For δ = 0.05 we obtain Therefore, we get γ = 1 and rc = 2.01 in the equilibration formula (7.2). Figure 3 shows three different samples of the diffusion coefficient and the corresponding FE approximations of the PDE solution. The FE discretization parameters are given by h = 0.4 ⋅ 2 −( −1) for l = 1, ..., 7. We set u ref = u 7,ε W ,ε l , where the approximation parameters ε W and ε l are choosen according to (7.2) and we compute M = 100 samples to estimate the strong error by the Monte Carlo estimator (see (7.3)). In this experiment we investigate the strong error convergence rate for the sample-adapted FE approach as well as convergence rate for the non-adapted FE approach (see Section 6). In Subsection 7.3.1 we described two approaches to simulate Poisson subordinators. We run this experiment with both approaches: first, we approximate the Poisson process via sampling on an equidistant level-dependent grid and, in a second run of the experiment, we simulate the Poisson subordinators exactly using the Uniform Method described in [30,Section 8.1.2]. The convergence results for the both approaches for this experiment are given in the Figure 4.
We see a convergence rate of approximately 0.7 for the standard FEM discretization and full order convergence (κ a ≈ 1) for the sample-adapted approach. On the right hand side of Figure  4 on sees that the sample-adapted approach is more efficient in terms of computational effort if we consider the error-to-(averaged)DOF-plot. Only on the first level the standard FEM approach seems to be more efficient (pre-asymptotic behaviour). If we compare the results for the approximation method with the Uniform Method (see 7.3.1), we find that, while the convergence rates are the same, the constant of the error in the sample-adapted approach is slightly smaller for the Uniform Method. This shift is exactly the additional error resulting from an approximation of the subordinators in the approximation approach. We also see that, compared to the approximation approach, on the lower levels the averaged degrees of freedom in the sample-adapted FEM approach is slightly higher if we simulate the Poisson subordinators exactly. This is caused by the fact that in this case we do not approximate the discontinuities of the field which are generated by the Poisson processes. This results in a higher average number of degrees on freedom on the lower levels because discontinuities are more likely close to each other.  Figure 5 shows samples of the diffusion coefficient and the corresponding FE approximation of the PDE solution. We estimate the strong error convergence rate for this problem in the same way as in the previous example using M = 250 samples and we use the approximation approach to simulate the Poisson subordinators (see Subsection 7.3.1). Convergence results are given in Figure 6. As in the experiment with mixed Dirichlet-Neumann boundary conditions we obtain convergence order of κ a ≈ 0.7 for the standard FEM approach and full order convergence for the sampleadapted approach. Also in case of homogeneous Dirichlet boundary conditions the sample-adapted FEM is more efficient in terms of the averaged number of degrees of freedom.
7.3.4. Poisson subordinators: high intensity and mixed boundary conditions. In this section we want to consider subordinators with higher intensity, resulting in a higher number of discontinuities in the diffusion coefficient. Therefore, we consider l 1 and l 2 to be Poisson (5) processes.
We set the cut-off level K of the subordinators in Equation (4.4) to K = 15. For this choice it is reasonable to expect that this cut-off has no numerical influence since for j = 1, 2. However, setting K = 15 means that we have to simulate the GRF W 2 on the domain [0, 15] 2 which would be time consuming. Therefore, we set K = 1 instead and consider the downscaled processesl for t ∈ [0, 1] and j = 1, 2. The variance parameter of the field W 1 is chosen to be σ 1 = 1 and the parameters of the GRF W 2 are set to be σ 2 = 0.3 and r 2 = 0.5. Figure 7 shows samples of the coefficient and the corresponding pathwise FEM solution. As in the first experiment, we again run this experiment using both methods described in Subsection 7.3.1: the approximation approach using Poisson-distributed increments and the Uniform Method. We use the discretization steps h = 0.1 ⋅ 1.7 −( −1) for = 1, . . . , 7 and M = 150 samples.
In Figure 8 we see that we get almost full order convergence for the sample-adapted FE method for both approximation approaches of the Poisson processes. Compared to the lowintensity examples with Poisson(1)-subordinators given in Subsection 7.3.2 and 7.3.3, we get a slightly lower convergence rate of approximately 0.55 for the standard FEM approach. This holds for both approximation methods of the Poisson subordinators. Hence, we see that the way how the Poisson-subordinators are simulated seems to have no effect on the convergence rate. 7.3.5. Poisson subordinators of a GRF with short correlation length: high intensity and mixed boundary conditions. In our construction of the jump-diffusion coefficient, the jumps are generated by the subordinated GRF. To be precise, the number of spatial jumps is determined by the subordinators and the jump intensities (in terms of the differences in height between the jumps) are essentially determined by the GRF W 2 . This fact allows to control the jump intensities of the diffusion coefficient by the correlation parameter of the underlying GRF W 2 . In the following experiment we want to investigate the influence of the jump intensities of the diffusion coefficient on the convergence rates. In Subsection 7.3.4 we subordinated a Matérn-1.5-GRF with pointwise standard deviation σ 2 = 0.3 and a correlation length of r 2 = 0.5. In the following experiment we set the standard deviation of the GRF W 2 to σ 2 = 0.5 and the correlation length to r 2 = 0.1 and leave all the other parameters unchanged. Figure 9 compares the resulting GRF with the field W 2 with parameters σ 2 = 0.3 and r 2 = 0.5 which we used in Subsection 7.3.4. Subordinating the GRF with small correlation length (right plots in Figure 9) result in higher jump intensities in the diffusion coefficient as the subordination of the GRF with higher correlation length (left plots in Figure 9). Figure 10 shows samples of the diffusion coefficient and the corresponding PDE solutions where the parameters of W 2 are σ 2 = 0.5 and r 2 = 0.1. As expected, the resulting jump coefficient shows jumps with a higher intensity compared to the jump coefficient in the previous experiment where we used the GRF W 2 with parameters σ 2 = 0.3 and r 2 = 0.5 (see Figure 7).
We estimate the strong error convergence rate using this high-intensity jump coefficient using M = 200 samples and approximate the Poisson subordinators by the Uniform Method. Figure 11 shows that for the GRF W 2 with small correlation length the convergence rates are reduced for both approaches: the standard FEM approach and the sample-adapted version. We cannot preserve full order convergence in the sample-adapted FEM but observe a convergence rate of approximately 0.75. In the non-adapted approach we obtain a convergence rate of approximately 0.45. Looking at the error-to-(averaged)DOF-plot on the right hand side of Figure 11 we see that still the sample-adapted approach is by a large margin more efficient in terms of computational effort. This experiment confirms our expectations since the FEM convergence rate has been shown to be strongly influenced by the regularity of the jump-diffusion coefficient (see e.g. [12] and [28]).

Gamma subordinators.
In order to also consider Lévy subordinators with infinite activity we take Gamma processes to subordinate the GRF in the remaining numerical examples. We set the standard deviation of the GRF W 1 to be σ 1 = 1.5 and we choose σ 2 = 0.3 and r 2 = 1 for the Matérn-1.5-GRF W 2 and leave the other parameters unchanged. For a G , b G > 0, a Gamma(a G , b G )-distributed random variable admits the density function where Γ(⋅) denotes the Gamma function. A Gamma process (X t , t ≥ 0) has independent Gamma distributed increments. Being precise, X t − X s ∼ Gamma(a G ⋅ (t − s), b G ) for 0 < s < t (see [30,Chapter 8]). The following lemma is essential to approximate the Gamma processes.
Proof. We prove the assertion by induction. The case n = 0 is trivial. For an arbitrary natural number n ≥ 0 we calculate where we used the Theorem of Bohr-Mollerup (see [6]) in the last step.
In our numerical experiments we choose l j to be a Gamma(4, 10)-process for j = 1, 2. Since increments of a Gamma process are Gamma-distributed random variables it is straightforward to generate values of a Gamma process on grid points ( We then use the piecewise constant extension of the simulated values {l j (x i ), i = 0, . . . , N l − 1, j = 1, 2} to approximate the Lévy subordinators: for j = 1, 2. Note that in this case Assumption 4.2 v is fulfilled with for any fixed η < +∞. To see that we consider a fixed s ∈ N with s ≤ η and calculate for an arbitrary x ∈ [0, 1) with x ∈ [x i , x i+1 ):  We set the diffusion cut-off to K = 2 since in this case we obtain for j = 1, 2. The use of infinite-activity Gamma subordinators in the diffusion coefficient does not allow anymore for a sample-adapted approach to solve the PDE problem. Hence, we only use the standard FEM approach to solve the PDE samplewise and estimate the strong error convergence. We use M = 200 samples to estimate the strong error on the levels = 1, . . . , 5 where we set the non-adaptive FEM solution u 7,ε W ,ε l on level L = 7 to be the reference solution. We choose the FEM discretization steps to be h = 0.4 ⋅ 2 −( −1) for = 1, . . . , 7. Figure 13 shows a convergence rate of approximately 0.8 for the standard-FEM approach. Since we do not treat the discontinuities in a special way we cannot expect full order convergence. In fact, the given convergence is comparably good since in general we cannot prove a higher convergence order than 0.5 for the standard deterministic FEM approach without special treatment of the discontinuities (see [7] and [12]). The convergence rate of approximately 0.8 in this example is based on the comparatively large correlation length of the underlying GRF W 2 (see 7.3.5).
In Subsection 7.3.5 we investigated the effect of a rougher diffusion coefficient on the convergence rate for Poisson(5)-subordinators. In the following experiment we follow a similar strategy and use a shorter correlation length in the GRF W 2 which is subordinated by Gamma processes. Therefore, we choose the parameters of the Matérn-1.5-GRF W 2 to be σ 2 = 0.3 and r 2 = 0.05. Figure 14 shows a comparison of the resulting GRFs W 2 with the different correlation lengths. In Figure 15, the GRF with small correlation length results in higher jumps of the diffusion coefficient and stronger deformations of the corresponding PDE solution compared to the previous example (see Figure 12). We estimate the strong error taking M = 200 samples where we use the non-adapted FEM solution u 9,ε W ,ε l on level L = 9 as reference solution and choose the FEM discretization steps to be h = 0.1 ⋅ 1.5 −( −1) for = 1, . . . , 9. Figure 16 shows the convergence on the levels = 1, . . . , 6.
We observe a convergence rate of approximately 0.45 which is significantly smaller than the rate of approximately 0.8 we obtained in the example where we used a GRF W 2 with correlation length r 2 = 1 (see Figure 13). This again confirms that, for subordinated GRFs, the convergence rate of the FE method is highly dependent on the correlation length of the underlying GRF W 2 Proof.
We estimate using Hölder's inequality: For the second factor we estimate using the independence of W 1 and W 2 P(ess sup x∈D a K (x) ≥ A) = 1 − P(ess sup x∈D a K (x) ≤ A) where we used Equations (A.4) and (A.7) in the fourth step and therefore we obtain