On properties and applications of Gaussian subordinated L\'evy fields

We consider Gaussian subordinated L\'evy fields (GSLFs) that arise by subordinating L\'evy processes with positive transformations of Gaussian random fields on some spatial domain $\mathcal{D}\subset \mathbb{R}^d$, $d\geq 1$. The resulting random fields are distributionally flexible and have in general discontinuous sample paths. Theoretical investigations of the random fields include pointwise distributions, possible approximations and their covariance function. As an application, a random elliptic PDE is considered, where the constructed random fields occur in the diffusion coefficient. Further, we present various numerical examples to illustrate our theoretical findings.

1. Introduction. Ample applications that are modelled stochastically require stochastic processes or random fields, which allow for discontinuities or possess higher distributional flexibility than the standard Gaussian model (see for example [27], [38] and [43]). In case of a one-dimensional parameter domain, where the parameter often represents time, Lévy processes are widely used (see e.g. [38] for some applications in finance). For higher-dimensional parameter domains, some extensions of the Gaussian model have been proposed in the literature: The authors of [19] consider (smoothed) Lévy noise fields with high distributional flexibility and continuous realizations in the context of random PDEs. In [10], the authors propose a general random field model which allows for spatial discontinuities with flexible jump geometries. However, in its general form, it not easy to investigate theoretical properties of the random field itself. In the context of Bayesian inversion, the level set approach combined with Gaussian random fields is often used as a discontinuous random field model (see, e.g. [28], [17] and [18]). The distributional flexibility of the resulting model is, however, again restricted since the stochasticity is governed by the Gaussian field. Another extension of the Gaussian model has been investigated in the recent paper [32]. The construction is motivated by the subordinated Brownian motion, which is a Brownian motion time-changed by a Lévy subordinator (i.e. a non-decreasing Lévy process). As a generalization, the authors consider Gaussian random fields on a higher-dimensional parameter domain subordinated by several independent Lévy subordinators. The constructed random fields allow for spatial jumps and display great distributional flexibility. However, the jump geometry is restricted in the sense that the jump interfaces are always rectangular.
In this paper, we investigate another specific class of discontinuous random fields: the Gaussian subordinated Lévy fields (GSLF). Motivated by the subordination of standard Lévy processes, the GSLF is constructed by subordination of a general real-valued Lévy process by a transformation of a Gaussian random field. It turns out that the resulting fields allow for spatial discontinuities with flexible jump geometries, are distributionally flexible and relatively easy to simulate. Besides a theoretical investigation of the constructed random fields, we present numerical examples and introduce a possible application of the fields in the diffusion coefficient of an elliptic model PDE. Such a problem arises, for instance, in models for subsurface flow in heterogeneous/porous media (see [10], [14], [34] and the references therein).
The rest of the paper is structured as follows: In Section 2 we shortly introduce Lévy processes and Gaussian random fields, which are crucial for the definition of the GSLF in Section 3. The pointwise distribution of the random fields are investigated in Section 4, where we derive a formula for its (pointwise) characteristic function. Section 5 deals with numerical approximations of the GSLF, including an investigation of the approximation error and the pointwise distribution of the approximated fields. Our theoretical investigations are concluded in Section 6, where we derive a formula for the covariance function. In Section 7, we present various numerical examples which  There exist parameters γ ∈ R, b ∈ R + and a measure ν on R such that the characteristic function φ l(t) of the Lévy process l admits the representation φ l(t) (ξ) ∶= E(exp(iξl(t))) = exp(tψ(ξ)), ξ ∈ R, for t ∈ T . Here, ψ denotes the characteristic exponent of l which is given by Further, the measure ν satisfies R min(y 2 , 1) ν(dy) < ∞, and is called Lévy measure and (γ, b, ν) Lévy triplet.

Gaussian random fields.
A random field defined over the Borel set D ⊂ R d is a family of real-valued random variables on the probability space (Ω, F, P) parametrized by x ∈ D. The Gaussian random field defines an important example. Definition 2.3. (see [3,Sc. 1.2]) A random field (W (x), x ∈ D) on a d-dimensional domain D ⊂ R d is said to be a Gaussian random field (GRF) if, for any x (1) , . . . , x (n) ∈ D with n ∈ N, the ndimensional random variable (W (x (1) ), . . . , W (x (n) )) follows a multivariate Gaussian distribution. In this case, we define the mean function µ W (x (i) ) ∶= E(W (x (i) )) and the covariance function q W (x (i) , x (j) ) ∶= E((W (x (i) ) − µ W (x (i) ))(W (x (j) ) − µ W (x (j) ))), for x (i) , x (j) ∈ D. The GRF W is called centered, if µ W (x) = 0 for all x ∈ D.
Note that every GRF is determined uniquely by its mean and covariance function. The GRFs considered in this paper are assumed to be mean-square continuous, which is a common assumption (cf. [3]). We denote by Q ∶ L 2 (D) → L 2 (D) the covariance operator of W which is defined by Q(ψ)(x) = D q W (x, y)ψ(y)dy, for x ∈ D, for ψ ∈ L 2 (D). Here, L 2 (D) denotes the Lebesgue space of all square integrable functions over D (see, for example, [2]). If D is compact, it is well known that 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 [42, Theorem VI.3.2 and Chapter II.3]). A GRF W is called stationary if the mean function µ is constant and the covariance function q W (x (1) , x (2) ) only depends on the difference x (1) − x (2) of the values x (1) , x (2) ∈ D and a stationary GRF W is called isotropic if the covariance function q W (x (1) , x (2) ) only depends on the Euclidean length x (1) − x (2) 2 of the difference of the values x (1) , x (2) ∈ D (see [3], p. 102 and p. 115).
Example 2.4. The Matérn-GRFs are a class of continuous GRFs which are commonly used in applications. For a certain smoothness parameter ν > 1 2, correlation parameter r > 0 and variance σ 2 > 0, the Matérn-ν covariance function is defined by Here, Γ(⋅) is the Gamma function and K ν (⋅) is the modified Bessel function of the second kind (see [22, Section 2.2 and Proposition 1]). A Matérn-ν GRF is a centered GRF with covariance function q M .
3. The Gaussian subordinated Lévy field. In this section we define Gaussian subordinated Lévy fields. Since their construction is motivated by the subordination of standard Lévy processes we shortly repeat this procedure: If l denotes a Lévy process and S denotes a Lévy subordinator (i.e. a non-decreasing Lévy process), which is independent of l, the time-changed process t ↦ l(S(t)), t ≥ 0, is called subordinated Lévy process. It can be shown that this process is again a Lévy process (cf. [6,Theorem 1.3.25]).
In order to construct the GSLF we consider a domain D ⊂ R d with 1 ≤ d ∈ N. Let l = (l(t), t ≥ 0) be a Lévy process, W ∶ Ω × D → R be an F ⊗ B(D)-measurable GRF which is independent of l and F ∶ R → R + be a measurable, non-negative function. The Gaussian subordinated Lévy field is defined by Note that assuming the GRF W to have continuous paths is sufficient to ensure joint measurability (see [5,Lemma 4.51]). Since the Lévy process l is in general discontinuous, the GSLF L has in general discontinuous paths. This is demonstrated in Figure 1, which shows samples of the GSLF. Remark 3.1. Discontinuous random fields often serve as prior model to solve inverse problems with a Bayesian approach as an alternative to the standard Gaussian prior (see e.g. [28,27]). In these situations, the prior model is often set to be a Gaussian(-related) level-set function. One way to construct such a prior model is as follow: where n ∈ N, (u i , i = 1, . . . , n) ⊂ R are fixed and with fixed levels (c i , i = 1, . . . , n) ⊂ R and a GRF W (see e.g. [18,17]). The GSLF, as defined above, may be interpreted as a generalization of the Gaussian level-set function.
4. The pointwise characteristic function of a GSLF. In the following section we derive a formula for the pointwise characteristic function of the GSLF, which determines the pointwise distribution entirely. Such a formula is especially valuable in applications, where distributions of a random field have to be fitted to data observed from real world phenomena. We start with a technical Lemma on the computation of expectations of functionals of the GSLF.
Lemma 4.1. Let l = (l(t), t ≥ 0) be a stochastic process with a.s. càdlàg paths and W + be a real-valued, non-negative random variable which is stochastically independent of l. Further, let g ∶ R → R be a continuous function. It holds Proof. The proof follows by the same arguments as in the proof of [8,Lemma 4.1].
Remark 4.2. Note that Lemma 4.1 also holds for complex-valued, continuous and bounded functions g ∶ R → C (cf. [8,Remark 4.2]). Further, we emphasize that Lemma 4.1 also holds for an R d + -valued random variable W + and a random field (l(t), t ∈ R d + ), independent of W + , which is a.s. càdlàg in each variable, i.e. for P-almost all ω ∈ Ω, it holds for t (n) j ↘ t j , for n → ∞, j = 1, . . . , d, and any t = (t 1 , . . . , t d ) ∈ R d + . With Lemma 4.1 at hand, we are able to derive a formula for the pointwise characteristic funciton of the GSLF. Corollary 4.3. Let l = (l(t), t ≥ 0) be a Lévy process with Lévy triplet (γ, b, ν) and W = (W (x), x ∈ D) be an independent GRF with pointwise mean µ W (x) = E(W (x)) and variance σ W (x) 2 ∶= V ar(W (x)) for x ∈ D. Further, let F ∶ R → R + be measurable. It holds for x ∈ D, where ψ denotes the characteristic exponent of l defined by Proof. We consider a fixed x ∈ D and use Lemma 4.1 with g(⋅) ∶= exp(iξ⋅) and W + ∶= F (W (x)) to calculate where m is defined through m(z) = E(exp(iξl(z)) = exp(z ψ(ξ)), z ∈ R + , by the Lévy-Khinchin formula (see Theorem 2.2). Hence, we obtain 5. Approximation of the fields. The GSLF may in general not be simulated exactly since in most situations it is not possible to draw exact samples of the corresponding GRF and the Lévy process. The question arises how the GSLF may be approximated and if the corresponding approximation error may be quantified. In this section we answer both questions. We prove an approximation result for the GSLF where we approximate the GRF and the Lévy process separately. To be more precise, we approximate the Lévy processes using a piecewise constant càdlàg approximation on a discrete grid (see e.g. [9] and the remainder of the current section). The GRF may be approximated by a truncated Karhuen-Loève-expansion or using values of the GRF simulated on a discrete grid, e.g. via Circulant Embedding (see, e.g. [10,1,39] resp. [23,24]). Naturally, we have to start with some assumptions on the regularity of the GRF and the approximability of the Lévy process. Assumptions of this type are natural and well known in different situations (see e.g. [10,33,13]). For simplicity, we consider centered GRFs in this subsection.
Assumption 5.1. Let W be a zero-mean GRF on the compact domain D. We denote by q W ∶ D × D → R the corresponding covariance function and by ((λ i , e i ), i ∈ N) the eigenpairs associated to the corresponding covariance operator Q as introduced in Section 2.2. In particular, (e i , i ∈ N) is an orthonormal basis of L 2 (D).
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 F ∶ R → R + is Lipschitz continuous and globally bounded by iii l is a Lévy process on [0, C F ] with Lévy triplet (γ, b, ν) which is independent of W . Further, we assume there exists a constant η > 1 and càdlàg approximations l (ε l ) of this process such that for every s ∈ [1, η) it holds with δ ∈ (0, 1] and a constant C l > 0 which may depend on s but is independent of t and ε l .
We continue with a remark on Assumption 5.1.
Remark 5.2. Assumtion 5.1 i is natural for GRFs and guarantees certain regularity properties for the paths of the GRF (see e.g. [10,33,13]). Equation (5.1) ensures that we can approximate the Lévy subordinators in an L s -sense. This can be achieved under appropriate assumptions on the tails of the distribution of the Lévy processes, see [9, Assumption 3.6, Assumption 3.7, Theorem 3.21] and [33,Section 7]). There are several results in the direction of condition (5.2). For example, in [30] and [15], the authors formulate general assumptions on the Lévy measure which guarantee Equation (5.2) and similar properties. Further, in [30] the authors explicitly derive the rate δ in Equation (5.2) for several Lévy processes. In [12,Proposition 2.3], an exact polynomial time-dependence of the absolute moments of a Lévy process under the assumption that the absolute moment of the Lévy process exists up to a certain order was proven. In order to illustrate (5.2), we present a short numerical example: for three different Lévy processes, we estimate E( l(t) s ) for t = 2 i with i ∈ {1, 0, −1, . . . , −16} using  We close this subsection with a remark on a possible way to construct approximations of Lévy processes.
Remark 5.3. One way to construct a càdlàg approximation l (ε l ) ≈ l is a piecewise constant extension of the values of the process on a discrete grid: assume . . , N l } are the values of the process l on this grid. We define Such a construction yields an admissible approximation in the sense of Assumption 5.1 iii for many Lévy processes (see e.g. [33,Section 7] for the case of Poisson and Gamma processes).
Note that the values of the process on the grid may be simulated easily for many Lévy processes using the independent and stationary increment property.

Approximation of the GRF.
In this subsection, we shortly introduce two possible approaches to approximate the GRF W. Both approaches are admissible for the construction of approximations of the GSLF in Subsection 5.2. In the following and for the rest of the section, we consider D = [0, D] d . Assumption 5.1 allows conclusions to be drawn on the spatial regularity of the GRF W . The next lemma on the approximation of the GRF W follows from the regularity properties of W (see [33,Lemma 4.4] and [13]).
be an approximation of the GRF W on the discrete grids G (ε W ) which is constructed by point evaluation of the random field W on the grid and multilinear interpolation between the grid points. Under Assumption 5.1 i, it holds for n ∈ [1, +∞): , where β and α are the parameters from Assumption 5.1 i.
It follows by the Karhunen-Loève-Theorem (see e.g. [3, Section 3.2]) that the GRF W admits the representation and the sum converges in the meansquare sense, uniformly in x ∈ D. The Karhunen-Loève expansion (KLE) motivates another approach to approximate the GRF W : For a fixed cut-off index N ∈ N, we define the approximation W N by Proof. It follows by [10, Theorem 3.8] that which finishes the proof.
In the following, we denote by W N an approximation of the GRF of W . The notation is clear if we approximate the GRF by the truncated KLE. If we approximate W by sampling on a discrete grid (cf. Lemma 5.4), then we use ε W = 1 N and abuse notation to write W N = W (1 N ) . Regardless of the choice between these two approaches, we thus obtain an approximation for N ∈ N and n ∈ [1, +∞), where R(N ) = N −γ if we approximate the GRF by discrete sampling and R(N ) = N − β 2 if we approximate it by the KLE approach.

Approximation of the GSLF.
We are now able to quantify the approximation error for the GSLF where both components, the Lévy process and the GRF, are approximated.
Theorem 5.6. Let Assumption 5.1 hold and assume an approximation W N ≈ W of the GRF is given, as introduced in Subsection 5.1. For a given real number 1 ≤ p < η and a globally Lipschitz for N ∈ N and ε l > 0, where δ is the positive constant from (5.2).
Proof. We calculate using Fubini's theorem and the triangle inequality Further, we use the Lipschitz continuity of g together with Lemma 4.1 to obtain for any x ∈ D for all x ∈ D and, hence, For the second summand we calculate using Lemma 4.1 and Remark 4.
For C F > t 1 ≥ t 2 ≥ 0, the stationarity of l together with Assumption 5.1 iii yields Overall, we obtain for any ω ∈ Ω the pathwise estimatẽ for any x ∈ D. We apply Hölder's inequality and Equation (5.4) to obtain for any x ∈ D. Finally, we use the Lipschitz contnuity of g to obtain: Overall, we end up with 5.3. The pointwise distribution of the approximated GSLF. In Section 4 we ivestigated the pointwise distribution of a GSLF and derived a formula for its pointwise characteristic function. In Section 5, we demonstrated how approximations of the Lévy process l and the underlying GRF W may be used to approximate the GSLF and quantified the approximation error. This is of great importance especially in applications, since it is in general not possible to simulate the GRF or the Lévy process on their continuous parameter domains. The question arises how such an approximation affects the pointwise distribution of the field. For this purpose, we prove in the following a formula for the pointwise charateristic function of the approximated GSLF.
Corollary 5.7. We consider a GSLF with Lévy process l and an independent GRF W . Let l (ε l ) ≈ l be a càdlàg approximation of the Lévy process and W N ≈ W be an approximation of the GRF, where we assume that the mean function µ W of the GRF W is known and does not need to be approximated. For x ∈ D, the pointwise characteristic function of the approximated GSLF is given by denotes the pointwise characteristic function of l (ε l ) . Further, consider the discrete grid If the Lévy process is approximated according to Remark 5.3 and the GRF is approximated with the truncated KLE, that is Here, σ 2 W,N is the variance function of the approximation W N , which is given by The function ψ denotes the characteristic exponent of l defined by and, for any positive real number x, we denote by ⌊x⌋ the largest integer smaller or equal than x.
Proof. As in the proof of Corollary 4.3, we use Lemma 4.1 to calculate In the next step, we compute the charateristic function of the approximation l (ε l ) of the Lévy process constructed as described in Remark 5.3. We use the independence and stationarity of the increments of the Lévy process to obtain where (l ε l k , k ∈ N) are i.i.d. random variables following the distribution of l(ε l ) and we denote by D = equivalence of the corresponding probability distributions. The convolution theorem (see e.g. [29,Lemma 15.11]) yields the representation for t ∈ [0, C F ) and ξ ∈ R. Therefore, we obtain The second formula for the characteristic function of the approximated field now follows from the 6. The covariance structure of GSLFs. In many modeling applications one aims to use a random field model that mimics a specific covariance structure which is, for example, obtained from empirical data. In such a situation it is useful to have access to the theoretical covariance function of the random fields used in the model. Therefore, we derive a formula for the covariance function of the GSLF in the following section.
Lemma 6.1. Assume W is a GRF and l is an independent Lévy process with existing first and second moment. We deonte by µ W , σ 2 W and q W , the mean, variance and covariance function of W and by µ l (t) = E(l(t)), µ (2) l (t) ∶= E(l(t) 2 ) the functions for the first second moment of the Lévy process l. For x ≠ x ′ ∈ D, the covariance function of the GSLF L is given by where we define . For x = x ′ ∈ D, the pointwise variance is given by Proof. We compute using Lemma 4.1, for x ∈ D µ L (x) ∶= E(L(x)) = E l(F (W (x))) = E µ l (F (W (x)))) .
Putting these results together, we end up with The assertion now follows from the fact that In this section, we present numerical experiments on the theoretical results presented in the previous sections. The aim is to investigate the results of existing numerical methods and to illustrate the theoretical properties of the GSLF which have been proven in the previous sections, e.g. the pointwise distribution of the approximated fields or the quality of this approximation (see Theorem 5.6). Further, the presented numerical experiments and methods may also be useful for fitting the GSLF to existing data in various applications.
7.1. Pointwise characteristic function. Corollary 4.3 gives access to the pointwise characteristic function of the GSLF L(x) = l(F (W (x))), x ∈ D, with a Lévy process l and a GRF W , which is independent of l. Using the characteristic function and the Fourier inversion (FI) method (see [21]) we may compute the pointwise density function of the GSLF. Note that in both of these steps, the application of Corollary 4.3 and of the Fourier inversion theorem, numerical integration is necessary which may be inaccurate or computationally expensive. In this subsection, we choose specific Lévy processes together with a GRF and compare the computed density function of the corresponding GSLF with the histogram of samples of the simulated field. To be precise, we choose a Matérn-1.5-GRF, a Gamma process, set F = ⋅ + 1 and consider the GSLF L(x) = l(F (W (x)) on D = [0, 1] 2 . The distributions and the corresponding density functions are presented in Figure  3. In line with our expectations, the FI approach perfectly matches the true distribution of the field at (1, 1). 7.2. Pointwise distribution of approximated fields. In Subsection 7.1 we presented a numerical example regarding the pointwise distribution of the GSLF. In applications, it is not possible to draw samples from the exact GSLF and, hence, one has to use approximations instead.  The effect of such an approximation on the pointwise distribution of the random field has been investigated theoretically in Subsection 5.3. In this section, we aim to provide a numerical example in order to visualize the distributional effect of sampling from an approximation of the GSLF. For a specific choice of the Lévy process l, the transformation function F and the GRF W , we use Corollary 5.7 and the FI method to compute the pointwise distribution of the approximated field l (ε l ) (F (W N (x))) ≈ l(F (W (x))) for different approximation parameters ε l (resp. N ) of the Lévy process (resp. the GRF). The computed densities are then compared with samples of the approximated field. We set D = [0, 1] 2 and consider the evaluation point x = (0.4, 0.6). The GRF W is given by the KLE where the eigenbasis is defined by with ν = 0.6 (see [20,Appendix A]). Further, we set F (x) = 1 + min( x , 30) and choose l to be a Gamma(3,10)-process. The threshold 30 in the definition of F is large enough to have no effect in our numerical example since the absolute value of W does not exceed 30 in all considered samples. We use piecewise constant approximations l (ε l ) of the Lévy process l on an equidistant grid with stepsize ε l > 0 (see Remark 5.3) and approximate the GRF W by the truncated KLE with varying truncation indices N . To be more precise, we choose 6 different approximation levels for these two approximation parameters, as described in Table 1. For each discretization level, we compute the characteristic function using Corollary 5.7, approximate the density function via FI and compare it with 10 5 samples of the approximated field evaluated at the point (0.4, 0.6). The results are given in Figure 4. As expected from Subsection 7.1, we see that the densities of the pointwise distribution of the approximated GSLF, which are approximated by FI, fit the actual samples of the approximated field accurately. In order to get a better impression of the influence of the approximation on the different levels, Figure 5 shows the densities of the evaluated approximated GSLF on the different levels, together with the density of the exact field. We see how the densities of the approximted GSLF converge to the density of the exact field for ε l → 0 and N → ∞. The results also show that the effect of the approximation of the GSLF should not be unterestimated: on the lower levels, we obtain comparatively large deviations of the pointwise densities from the density of the exact field, which should be taken into account in applications. Obviously, the effect of the approximation depends heavily on the specific choice of the Lévy process and the underlying GRF.
7.3. Numerical approximation of the GSLF. In Section 5, we considered approximations l (ε l ) ≈ l of the Lévy process and W N ≈ W of the GRF and derived an error bound on the corresponding approximation l (ε l ) (F (W N (x))) ≈ l(F (W (x))) (see Theorem 5.6). In fact, under Assumption 5.1, if we choose the approximation parameter N of the GRF W such that R(N ) ∼ ε In this section, we present numerical experiments to showcase this approximation result.  λ k e k (x, y)Z k , be a GRF with corresponding eigenbasis With this choice, Assumption 5.1 i is satisfied with α = 1 and 0 < β < 2ν − 1 , since where we used that λ k ≤ Ck −2ν for k ∈ N. Hence, we obtain by Lemma 5.5 for 0 < β < 2ν − 1, i.e. R(N ) = N − β 2 in Equation (5.4) and Theorem 5.6. In our experiments, we choose the Lévy subordinator l to be a Poisson or Gamma process. Approximations l (ε l ) ≈ l of these processes satisfying of Assumption 5.1 iii may be obtained by piecewise constant extensions of values of the process on a grid with stepsize ε l (see Remark 5.3). For these processes, we obtain δ = 1 in (5.2) (see Remark 5.2 and [33, Section 7]). Overall, Theorem 5.6 yields the error bound for 0 < β < 2ν − 1. Hence, in order to equilibrate the error contributions from the GRF approximation and the approximation of the Lévy process, one should choose the cut-off index N of the KLE approximation of the GRF W according to In our experiments, we set the approximation parameters of the Lévy process to be ε l = 2 − , for = 2, . . . , 10 and the cut-off indices of the GRF are choosen according to (7.2). In order to verify (7.3), we draw 150 samples of the random variable l (ε l ) (F (W N )) − l(F (W )), for the described approximation parameters to estimate the corresponding L p norm of the approximation error for p ∈ {1, 2, 2.5, 3, 3.3} which is expected to behave as O(ε 1 p l ). In our first example, we set l to be a Poisson(8)-process and ν = 2.5 in Equation (7.1). For each sample, we compute a reference field by taking 150 summands in the KLE of the GRF W and use an exact sample of a Poisson process computed by the Uniform Method (see [38,Section 8.1.2]). The resulting estimates for the L p approximation error are plotted against ε l for p ∈ {1, 2, 2.5, 3.3}, which is illustrated in the Figure 6. As expected, the approximated L p errors converge asymptotically with rate ε 1 p l for the considered moments p. The bottom plot of Figure 6 shows an overview of the L p approximation errors for different values of p and illustrates that the L p -errors indeed converge with different rates.
In our second numerical example, we set l to be a Gamma(2,4)-process. Further, we set ν = 2 in Equation (7.1). The approximations l (ε l ) ≈ l of the Lévy process on the different levels are again computed by piecewise constant extensions of values of the process on an equidistant grid with stepsize ε l . Aiming to verify Equation (7.3), we use 150 samples to estimate the L p approximation error, for p ∈ {1, 2, 2.5, 3, 3.3}. In order to compute a sufficiently accurate reference field in each sample, we use 500 summands in the KLE approximation of the GRF W and the Gamma process is computed on a reference grid with stepsize ε l = 2 −13 . The convergence of the estimated L p error plotted against the discretization parameter ε l is visualized in Figure 7, which shows the expected behaviour of the estimated L p error for the considered moments p. As in the previous experiment, we provide a plot with all estimated L p errors in one figure (see bottom plot in Figure  7), which again confirms that the L p -error of the approximation converge in ε l with rate 1 p for the considered values of p. 7.4. Empirical estimation of the covariance of the GSLF. Lemma 6.1 gives access to the covariance function of the GSLF. This is of interest in applications, since it is often required that a random field mimics a specific covariance structure which is determined by (real world) data. In this example, we choose specific GSLFs and spatial points to compute the corresponding covariance using Lemma 6.1 and compare it with empirically estimated covariances using samples of the GSLF. We choose D = [0, 2] 2 and estimate the covariance with a single level Monte Carlo estimator using a growing number of samples of the field L evaluated at the two considered points x, x ′ ∈ D. If M denotes the number of samples used in this estimation and E M (x, x ′ ) denotes the Monte Carlo estimation of the covariance, the convergence rate of the corresponding RMSE for a growing number of samples is −1 2, i.e.
In our experiments, we use Poisson and Gamma processes and W is set to be a centered GRF with squared exponential covariance function, i.e.
with pointwise variance σ 2 > 0, correlation length r > 0 and F = ⋅ . In our first experiment we choose l to be a Gamma process with varying parameters. The RMSE is estimated for a growing number of samples M and 100 independent MC runs. The results are shown in Figure 8. As expected, the approximated RMSE convergences with order O(1 √ M ) for both experiments. In our second example we set l to be a Poisson process. The results are shown in Figure 9. As in the previous example, we see a convergence rate of order O(1 √ M ) of the MC estimator for the covariance in each experiment. For small sample numbers, the error values for the estimation of the covariances might seem quite high: for example, in the left plot of Figure 9, using M = 100 samples, we obtain an approximated RMSE which is approximately two times the exact value of the covariance. This seems to large at first glance. However, one has to keep in mind that the RMSE is bounded by std(L(x) ⋅ L(x ′ )) √ M and the standard deviation of the product of the evaluated field may become quite large. In fact, in the left hand side of Figure 9 we obtain std(L(x) ⋅ L(x ′ )) ≈ 138, 46 which fits well to the observed results, being approximately 10 = √ 100 times the approximated RMSE for M = 100 samples. In practice, it is therefore important to keep in mind that the estimation of the covariance based on existing data might require a large number of observations. 8. GSLFs in elliptic PDEs. In the previous sections, we considered theoretical properties of the GSLF and presented numerical examples to validate and visualize these results. In the last section of this paper we present an application of the GSLF in the context of PDEs. This might be interesting, for example, to model subsurface flow in uncertain heterogeneous or fractured media. In such a modeling situation the media is often modeled by a random field, which should therefore be distributionally flexible and allow for spatial discontinuities, both of which are fulfilled in the case of the GSLF. This motivates the investigation of a general elliptic PDE where the GSLF occurs in the diffusion coefficient. We introduce the considered elliptic PDE in Subsection 8.1 following [31] and [33]. Spatial discretization methods are discussed in Subsection 8.2 and we conclude this section with numerical experiments in Subsection 8.3.
• W 1 and W 2 are zero-mean GRFs on D with P − a.s. continuous paths.
• l is a Lévy process on [0, C F ]. Note that we consider the case of a homogeneous Dirichlet boundary condition on Γ 1 in the theoretical analysis to simplify notation. Non-homogeneous Dirichlet boundary conditions could also be considered, since such a problem can always be formulated as a version of (8.1) -(8.3) with modified source term and Neumann data (see also [10, Remark 2.1]).
Next, we shortly introduce the pathwise weak solution of problem (8.1) -(8.3) following [33] and [31]. We denote by H 1 (D) the Sobolev space and by T the trace operator [16]). Further, we introduce the solution space where we take over the standard Sobolev norm, i.e. ⋅ V ∶= ⋅ H 1 (D) . We identify H with its dual space H ′ and work on the Gelfand triplet V ⊂ H ≃ H ′ ⊂ V ′ . Multiplying the left hand side of Equation (8.1) by a test function v ∈ V and integrating by parts (see e.g. [40, Section 6.3]) we obtain the following pathwise weak formulation of the problem: For fixed ω ∈ Ω and given mappings f (ω, ⋅) ∈ V ′ and g(ω, for all v ∈ V . The function u(ω, ⋅) is then called pathwise weak solution to problem (8.1) -(8.3) and the bilinear form B a(ω) and the linear operator F ω are defined by where the integrals in F ω are understood as the duality pairings: Remark 8.1. Note that any real-valued, càdlàg function f ∶ [a, b] → R on a compact interval [a, b], with a < b, is bounded. Otherwise, one could find a sequence (x n , n ∈ N) ⊂ [a, b] such that f (x n ) > n for all n ∈ N. In this case, since [a, b] is compact, there exists a subsequence (x n k , k ∈ N) ⊂ [a, b] with a limit in [a, b], i.e. x n k → x * ∈ [a, b] for k → ∞. Since f is rightcontinuous in x * and lim x↗x * f (x) exists, f is bounded in a neighborhood of x * , i.e. there exists δ, C > 0 such that f (x) ≤ C for x ∈ [a, b] with x − x * < δ, which contradicts the fact that f (x n k ) > n k since x n k − x * < δ and n k > C both hold for k large enough.
8.2. Spatial discretization of the elliptic PDE. In the following, we briefly describe numerical methods to approximate the pathwise solution to the random elliptic PDE, partially following [33, Section 6]. Our goal is to approximate the weak solution u to problem (8.1) -(8.3) with diffusion coefficient a given by Equation (8.4). Therefore, for almost all ω ∈ Ω, we aim to approximate the function u(ω, ⋅) ∈ V such that for every v ∈ V . We compute an approximation of the solution using a standard Galerkin approach with linear basis functions. Therefore, assume a sequence of finite-element subspaces is given, which is denoted by V = (V , ∈ N 0 ), where V ⊂ V are subspaces with growing dimension dim(V ) = d . We denote by (h , ∈ N 0 ) the corresponding refinement sizes with h → 0, for → ∞. Let ∈ N 0 be fixed and denote by {v where the coefficient vector c = (c 1 , . . . , c d ) T ∈ R d is determined by the linear system of equations ∈ N 0 , we denote by h ∶= max K∈K diam(K) the maximum diameter of the triangulation K and 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 there exists a positive regularity parameter κ a > 0 such that for P−almost all ω ∈ Ω it holds u(ω, ⋅) ∈ H 1+κa (D) and u L 2 (Ω;H 1+κa (D)) ≤ C u for some constant C u > ∞, we immediately obtain the following bound by Céa's lemma (see [10,Section 4], [33,Section 6], [26,Chapter 8]) u − u L 2 (Ω;V ) ≤ C h min(κa,1) , (8.7) for some constant C which does not depend on h . For general deterministic elliptic interface problems, one obtains a discretization error of order κ a < 1 and, hence, one cannot expect the full order of convergence κ a = 1 in general for standard triangulations of the domain (see [7] and [10]). Therefore, we present one possible approach to enhance the performance of the FE method for the considered elliptic PDE in Subsection 8.2.2. We point out that it is not possible to derive optimal rates κ a > 0 such that (8.7) holds for our general random diffusion coefficient (see also [33], [10]). However, we investigate the existence of such a constant numerically in Section 8.3. We close this subsection with a remark on the practical simulation of the GRFs W 1 , W 2 and the Lévy process l in the diffusion coefficient (8.4).
Remark 8.3. It is in general not possible to draw exact samples of the Lévy process and the GRFs in the diffusion coefficient (8.4). In practice, one has to use appropriate approximations l (ε l ) ≈ l of the Lévy process and W N 1 ≈ W 1 , W N 2 ≈ W 2 of the GRFs with approximation parameters ε l > 0 and N ∈ N, as introduced in Section 5. In the context of FE approximations of the PDE (8.1) -(8.3) with diffusion coefficient (8.4), the question arises, how to choose these approximation parameters in practice, given a specific choice of the FE approximation parameter h in (8.7). Obviously, the choice of ε l and N should depend on the FE parameter h , since a higher resolution of the FE approximation will be worthless if the approximation of the diffusion coefficient is poor.
In [33], the authors considered the PDE (8.1) -(8.3) with a different diffusion coefficient a. The coefficient considered in the mentioned paper also incorporates GRFs and Lévy processes, which in turn have to be approximated in practice, leading to an approximationã ≈ a. The authors derived a rule on the choice of the approximation parameters, such that the error contributions from the approximation of the diffusion coefficient and the FE discretization are equilibrated (see [33,Section 7.1]). We want to emphasize that this result is essentially based on the quantification of the approximation errorã − a of the corresponding diffusion coefficient in an L p (Ω, L p (D))-norm, for some p ≥ 1 (see [33,Theorem 4.8 and Theorem 5.11]).
In Theorem 5.6, we derived an error bound on the approximation error of the GSLF in the L p (Ω, L p (D))-norm under Assumption 5.1. A corresponding error bound on the approximation of the diffusion coefficient defined in (8.4) immediately follows under mild assumptions on Φ 1 and Φ 2 . Therefore, following exactly the same arguments as in [33] together with Theorem 5.6, we obtain the following rule for the practical choice of the approximation parameters ε l and N such that the overall error contributions from the approximation of the GRFs, the Lévy process and the FE discretization are equilibrated and the error is dominated by the FE refinement parameter h : For ∈ N, choose ε l and N such that For example, if we approximate the GRF by the KLE approach (see Subsection 5.1), one should choose the cut-off index such that N ≃ h − 4κa βδ with β from Assumption 5.1 i.

Adaptive finite elements.
As we pointed out in the last section, one cannot expect full order convergence (κ a = 1) for the FE method with linear elements in the considered elliptic problem due to the discontinuities in the diffusion coefficient. One common way to improve the FE method is to use triangulations which are adapted to the jump discontinuities in the sense that the spatial jumps lie on edges of elements of the triangulation. This approach leads to sampledependent triangulations which has been proven to enhance the performance of the FE method significantly compared to the use of standard triangulations, which are not adjusted to the jumps (see for example [10] and [33] and the references therein). In the cited papers, the jump locations are known and the jump geometries allow for an (almost) exact alignment of the triangulation to the interfaces due to their specific geometrical properties. This is, however, not the case for the diffusion coefficient considered in the current paper, where the spatial jump positions are not known explicitly, nor is it possible to align the triangulation exactly in the described sense due to the complex jump geometries.
Another possible approach to improve the FE method are adaptive Finite Elements (see e.g. [4] and [25] and the references therein). The idea is to identify triangles with a high contribution to the overall error of the FE approximation, which are then refined. For a given FE approximation u ∈ V ⊂ V of the solution u ∈ V , the discretization error u − u is estimated in terms of the approximated solution u and no information about the true solution u is needed. This may be achieved by the use of the Galerkin orthogonality and partial integration. For a given triangulation T = {K, K ⊂ D} with corresponding FE approximation u , the approximation error of the FE approximation to the considered elliptic problem (8.1) -(8.3) may be estimated by err T (u , a, f, g) = K∈T err T (u , a, f, g, K) 1 2 .
Here, err T (u , a, f, g, K) corresponds the error contribution of the approximation on the element K ∈ T , which may be approximated by where we denote by e ∈ ∂K an edge of the triangle K ∈ T and J(K, e) ∈ T is the unique triangle which is the direct neighbor of the triangle K along the edge e. Further, h K is the diameter of the triangle K and h e is the length of the edge e (see [25] and [11]). The triangle-based estimated error obtained by Equation (8.8) allows for an identification of the triangles with the largest contribution to the overall error. This may then be used to perform a local mesh refinement, which usually consists in the refinement of triangles with high error contribution. One common strategy is to start with an initial triangulation, compute the error contribution of each triangle according to Equation (8.8) and refine all triangles which have an approximated error contribution which is at least 50% of the error contribution with the largest approximated error (see, e.g., [41,Section 5]). where l is a Poisson(2)-process and W 1 is a centered GRF with squared exponential covariance function. In order to illustrate the adaptive mesh generation decribed above, we consider 3 samples of the diffusion coefficient and compute the adaptive meshes, where we use an initial mesh with FE refinement parameter h = 0.025 and refine all triangles which have an estimated error contribution exceeding 50% of the largest estimated error. Figure 10 shows the three samples of the diffusion coefficient and the corresponding adaptive triangulations. It is nicely illustrated, how the elementwise a-posteriori estimation of the error according to Equation (8.8) enables an identification of the triangles which lie near the discontinuities of the diffusion coefficient and, hence, have a comparatively high error contribution. This results in a local refinement leading to a higher mesh resolution near the discontinuity.

8.3.
Numerical experiments for the random elliptic PDE. In this section, we present numerical examples for the experimental investigation of the existence of a positive parameter κ a > 0 such that Equation (8.7) holds. Further, in the presented examples, we compare the performance of a standard FE discretization with the FE approach using the local refinement strategy described in Subsection 8.2.2. We consider the domain D = (0, 1) 2 and set f ≡ 10, a ≡ 0.1, Φ 1 = 0.1 exp(⋅), Φ 2 = ⋅ , F = min( ⋅ , 30) in Equation (8.4). W 1 and W 2 are centered GRFs with squared exponential covariance function (see Subsection 7.4), where we set σ = 0.5 and r = 1. We use the circulant embedding method (see cf. [23] and [24]) to draw samples of the GRFs W 1 , W 2 on an equidistant grid with stepsize h W = 1 200 = 0.005 and obtain evaluations everywhere on the domain by multilinear interpolation. Due to the high spatial regularity of W 1 and W 2 and the high correlation length, this stepsize is fine enough to ensure that the approximation error of the GRFs W 1 , W 2 is negligible. The Lévy process l is set to be a Poisson process with intensity parameter λ > 0, which will be specified later. It follows that l(t) − l(s) ∼ P oiss(λ(t − s)) for 0 ≤ s ≤ t and we draw samples from the Poisson process by the Uniform Method (see [   We use a reference grid on D with 801 grid points for interpolation and prolongation and take u ref ∶= u 8 as the pathwise reference solution. The RMSE is estimated for the standard FE method and for the FE method with adaptive refinement as described in Subsection 8.2.2. In order to obtain comparable approximations on each FE level, we compute the adaptive meshes as follows: for ∈ {1, . . . , 5}, we denote by n the number of triangles in the non-adaptive triangulation with FE mesh refinement parameter h . The (sample-dependent) adaptive mesh on level is obtained by performing the local refinement procedure described in Subsection 8.2.2 until the number of triangles in the adaptive mesh exceeds n . The resulting mesh is then used to compute the adaptive FE approximation on level .
8.3.1. Homogeneous Dirichlet boundary conditions. In our first experiment, we choose homogeneous Dirichlet boundary contitions on ∂D and set λ = 2. Figure 11 shows samples of the diffusion coefficient and the corresponding PDE solutions approximated by the FE method. We use M = 150 samples so approximate the RMSE according to Equation (8.9) with the non-adaptive and the adaptive FE approach. Figure 12 shows the approximated RMSE plotted against the inverse FE refinement parameter. For the adaptive FE approximations, the sample-dependent mesh refinement parameters on each discretization level are averaged over all 150 samples. We see a convergence with rate ≈ 0.65 for the standard FE method, which is in line with our expectations (see Section 8.2). Further, we observe that the adaptive refinement strategy yields an improved convergence rate of ≈ 0.85 and a smaller estimated RMSE on the considered levels. The right hand side of Figure 12 shows the estimated RMSE plotted against the degrees of freedom of the linear FE-system on each discretization level. For the adaptive FE method, the degrees of freedom are averaged over all samples. After a pre-asymptotic behaviour on the first and second level, we see that the adaptive FE method performs more efficient in terms of the degrees of freedom: the adaptive FE method achieves a certain RMSE with less degrees of freedom compared to the standard FE method. for ω ∈ Ω. Further, we set λ = 3, which results in a higher number of jumps in the diffusion coefficient. Figure 13 shows samples of the diffusion coefficient and the corresponding PDE solutions approximated by the FE method. As in the first experiment, we use M = 150 samples to approximate the RMSE according to (8.7) with the standard FE approach and the adaptive approach. The approximated values are plotted against the inverse FE refinement parameter. The results are presented in Figure 14, which shows a convergence rate of ≈ 0.6 for the standard FE approach, which is slightly smaller than the observed convergence rate of ≈ 0.65 in the first numerical example, where we imposed homogeneous Dirichlet boundary conditions and considered a Poisson process with intensity parameter λ = 2, leading to a smaller expected number of jumps in the diffusion coefficient. Further, Figure 14 shows a convergence rate of ≈ 0.85 for the adaptive FE approach and smaller magnitudes of the RMSE compared to the standard FE method. The right plot of Figure 14 reveals the higher efficiency of the adaptive FE method compared to the standard approach in the sense that the number of degrees of freedom, which are necessary to achieve a certain error, is significantly smaller compared to the standard FE approach. Further, we see that the performance difference of the standard and the adaptive FE method is larger compared to the first example due to the higher number of expected jumps in the diffusion coefficient, which renders the adaptive local refinement strategy even more suitable for this problem. Overall, the results are in line with our expectations.

Mixed
Dirichlet-Neumann boundary conditions and jump-accentuated coefficients. In our last experiment we consider mixed Dirichlet-Neumann boundary conditions as in the previous section. The diffusion coefficient is set to be a(x, y) = 0.01 + 0.01 exp(W 1 (x, y)) + 30 l(F (W 2 (x, y)), for (x, y) ∈ D, where l is a Poisson(4)-process. This leads to a jump-accentuated coefficient with high contrast. We take M = 500 samples to estimate the RMSE for the standard FE method and the adaptive approach according to (8.9). The results are presented in Figure 15. We obtain a convergence rate of ≈ 0.55 for the standard FE approach and ≈ 0.8 for the adaptive FE method, which is slightly slower than the observed rates in the previous experiment. This is in line with our expectations since we have an increased jump intensity in the underlying Poisson process and larger jump heights in the diffusion coefficient, both having a negative influence on the convergence rate of the FE method (see also [33,10,35]). Figure 15 also reveals, as expected, that the magnitude of the RMSE is significantly smaller in the adaptive FE approach. The right plot of Figure  15 demonstrates that the adaptive refinement strategy is able to achieve a certain error with significantly less degrees of freedom in the corresponding linear FE system, than the standard FE approach.