Elena Di Bernardino
email: elenadb@unice.fr

Anne Estrade
email: anne.estrade@u-paris.fr

Thomas Opitz
email: thomas.opitz@inrae.fr

Spatial extremes and stochastic geometry for Gaussian-based peaks-over-threshold processes

Keywords: Excursion sets, Extremal coefficient, Gaussian Kinematic Formula, Lipschitz-Killing Curvatures, Pareto Process, Regular variation

come

Introduction

Extremes of stochastic processes defined over Euclidean space R N for integer dimension N ≥ 1 (also called random fields in the following) have been extensively studied from the theoretical and statistical perspective. The approach of stochastic geometry of random fields focuses on geometric features of such fields. The study of the excursion sets of the process above some threshold level takes an important place since their geometric properties, such as the area, perimeter and number of connected components, provide relevant summaries of the spatial structure of the process. They are particularly relevant to gain better understanding of the clustering structure and dependence among extreme values by considering high thresholds. In this work, we will focus on the exceedance region within a nonempty domain S ⊂ R N of the stochastic process X = {X(s) : s ∈ S}. We make the assumption that X is almost surely of class C 3 , i.e., the 1 paths of the process are almost surely three times continuously differentiable. Given a real number u, the exceedance region of X in S above the threshold u (also called excursion set), is given by A(u, X, S) = {s ∈ S : X(s) ≥ u}.

(1)

Many results on the geometrical features of the random set in [START_REF] Abaach | Testing marginal symmetry of digital noise images through the perimeter of excursion sets[END_REF] have been established in the case of Gaussian processes X, such as the distribution of the maximum over space or the Euler-Poincaré characteristic describing the topological structure of excursion sets (see, e.g., the monographs of [START_REF] Robert | The geometry of random fields[END_REF] for an early account and of [START_REF] Robert | Random fields and geometry[END_REF] for a recent overview). Gaussian processes do not arise naturally as limits for extreme values. In extreme-value theory, pointprocess limits with each point defining a random function on R N take a pivotal role, and the well-studied max-stable processes arise from such point processes by taking locationwise maxima over all points. The importance of these processes is due to the fact that they are the only possible limits for linearly rescaled identically distributed (i.i.d.) copies X i , i = 1, 2, . . . , n of a stochastic process X, where the linearly rescaled processes (X i -b n )/a n are obtained using two functions a n > 0 and b n . A constructive characterization of the limit point processes with atoms Y i , i ∈ N, is given by the so-called spectral representation ( [START_REF] De | A spectral representation for max-stable processes[END_REF][START_REF] Schlather | Models for stationary max-stable random fields[END_REF][START_REF] De Haan | Extreme value theory: an introduction[END_REF])

{Y i (s), s ∈ R N , i ∈ N} = {Λ i V i (s), s ∈ R N , i ∈ N}, (2)
with the scaling variables Λ i = 1/U i for 0 < U 1 < U 2 < . . . the points of a unit-rate Poisson process on [0, ∞), and with random functions V i that are i.i.d. copies of some random function V on R N satisfying the moment constraint E[V + (s)] = 1, where a + = max(a, 0), for a ∈ R. Moreover, the series of processes {V i } must be independent of {Λ i }. Then, the sequence of processes {Y i } forms a Poisson point process (with points being functions) and consequently, the same holds for the sequence of variables {Y i (s)} for any fixed s (with points being scalars). The most widely used statistical extreme-value models for spatial processes (usually for geographic space R 2 ) are based on choosing V according to certain transformed Gaussian processes ( [START_REF] Anthony | Statistics of extremes[END_REF][START_REF] De | High-dimensional peaks-over-threshold inference[END_REF]). A first possibility is to choose V as the power of a centered Gaussian process ( [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF]), yielding the maxstable extremal-t processes. A second possibility consists of choosing V as a log-Gaussian process ( [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]), yielding a class of max-stable processes that extend the Brown-Resnick processes, for which V is log fractional Brownian motion [START_REF] Bruce | Extreme values of independent stochastic processes[END_REF]; we refer to this class as processes of Brown-Resnick type. These two model classes arise as limits when choosing X as certain scale or location mixtures of a Gaussian process: in a Gaussian location mixture (resp. scale mixture), the mean (resp. the standard deviation) of the Gaussian process is a random variable. Equation [START_REF] Adler | Excursion sets of three classes of stable random fields[END_REF] shows that the mixture structure persists in the construction of the extreme-value limit. Specifically, to obtain the location mixture representation, we consider log Y i (s) = log Λ i + log V i (s) when using a log-Gaussian process V . Instead of considering componentwise maxima, we can also study the process X conditional to an exceedance of its maximum taken over a compact domain S. Then, the existence of the point-process limit (and of the max-stable limit for maxima) is equivalent to the existence of a so-called Pareto process arising as the limit for the rescaled conditional process u -1 X | max s∈S X(s) > u as the threshold u tends to infinity ( [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF]). Specifically, Pareto process limits arise for Gaussian mixtures when the scale variable is regularly varying ( [START_REF] Huser | Bridging asymptotic independence and dependence in spatial extremes using Gaussian scale mixtures[END_REF]), or when the exponential of the location variable is regularly varying ( [START_REF] Krupskii | Factor copula models for replicated spatial data[END_REF]).

In general, Gaussian mixture processes possess stronger extremal dependence than Gaussian processes since joint exceedances above high levels occur more frequently. In contrast to asymptotic independence in Gaussian processes ( [START_REF] Sibuya | Bivariate extreme statistics[END_REF]), the mixture processes can be asymptotically dependent, meaning that the conditional probability P(F (X(s 1 )) > u | F (X(s 2 )) > u), with F the cumulative distribution function of X(s 1 ) and X(s 2 ) for s 1 = s 2 , has a positive limit as the threshold u tends to one. Extreme-value limit processes are characterized by asymptotic stability properties, and they possess asymptotic dependence when using Gaussian or log-Gaussian processes in the max-stable construction, except for the case of perfect negative Gaussian correlation. This means that a Pareto process conditioned on a threshold exceedance of its maximum has the same distribution as the original Pareto process up to rescaling. For cases where these asymptotic stability properties are too restrictive in practice, more general location or scale mixture constructions have recently been proposed to allow for more flexible statistical modeling of multivariate and spatial extremes ( [START_REF] Opitz | Modeling asymptotically independent spatial extremes based on Laplace random fields[END_REF], [START_REF] Huser | Bridging asymptotic independence and dependence in spatial extremes using Gaussian scale mixtures[END_REF], [START_REF] Krupskii | Factor copula models for replicated spatial data[END_REF], [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF]).

In spatial extreme-value analysis ( [START_REF] Anthony C Davison | Statistical modeling of spatial extremes[END_REF]), approaches based on threshold exceedances have been extensively studied recently ( [START_REF] Huser | Advances in statistical modeling of spatial extremes[END_REF]). However, summary statistics used in practice are typically defined from bivariate or low-dimensional observation vectors with respect to spatial distance between the locations, such as the bivariate extremal coefficient function or the tail correlation function ( [START_REF] Anthony C Davison | Geostatistics of dependent and asymptotically independent extremes[END_REF][START_REF] Strokorb | Tail correlation functions of max-stable processes[END_REF]). This restriction makes sense for the use with observations over irregularly spaced measurement locations (such as fixed weather stations on the ground). However, over the past decades the availability of datasets defined over regular and relatively dense spatial grids has strongly increased, and such data can provide an almost continuous cover of space. Data provided for a grid support arise from remote sensing techniques and from the output of many physical models (e.g., models for climate processes, fluid mechanics or population dynamics). In this context of spatial extreme-value analysis of gridded datasets, where often many temporal replicates are available, it makes sense to more systematically develop and use thresholdbased summary statistics of higher-than-bivariate order, such as by considering properties of excursion sets. This approach allows characterizing the joint behaviour over the full spatial domain and conveys a more complete picture of the extreme-value properties.

Here, we will take advantage of results on the stochastic geometry of Gaussian excursion sets and generalize them to provide novel results for more general Gaussian mixture processes, and for the corresponding Pareto limit processes with dependence of extremal-t or Brown-Resnick type. Due to peaks-over-threshold stability, the geometric properties of excursions above level u will not depend on u.

For geometric summaries of the excursion sets in (1) for N -dimensional Euclidean space, we focus on the so-called Lipschitz-Killing curvatures (LKCs) L k , k = 0, 1, . . . , N . Loosely speaking, L 0 is the Euler-Poincaré characteristic (i.e., the number of connected components for N = 1, or the difference between the number of connected components and the number of holes for N = 2), L N -1 is half the (N -1)-dimensional volume of the boundary (e.g., the perimeter for N = 2, the surface area for N = 3), and L N is the N -dimensional volume (i.e., the area for N = 2).

The expected values of LKCs of the excursion set have been studied in a wide variety of contexts (see, e.g., [START_REF] Robert | Random fields and geometry[END_REF][START_REF] Adler | Topological complexity of smooth random functions[END_REF][START_REF] Adler | Rotation and scale space random fields and the Gaussian kinematic formula[END_REF] for a focus on the Gaussian kinematic formula and [START_REF] Biermé | On the perimeter of excursion sets of shot noise random fields[END_REF] in the case where X is a shotnoise field). The expected value of L 0 was studied by [START_REF] Adler | Excursion sets of three classes of stable random fields[END_REF] for subGaussian random fields, a special case of Gaussian scale mixture arising for scale variables defined from stable distributions. Estimators for LKCs have also been studied and specific asymptotic results have been established. In the Gaussian framework, [START_REF] Estrade | A Central Limit Theorem for the Euler characteristic of a Gaussian excursion set[END_REF][START_REF] Di Bernardino | A test of Gaussianity based on the Euler characteristic of excursion sets[END_REF] studied the Euler-Poincaré characteristic, whereas the area (also called the sojourn time) is studied in [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF][START_REF] Pham | On the rate of convergence for central limit theorems of sojourn times of Gaussian fields[END_REF]. Central-limit theorems for LKCs were proposed by [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] and [START_REF] Müller | A central limit theorem for Lipschitz-Killing curvatures of Gaussian excursions[END_REF]. The setting of a stationary isotropic Gaussian field on R 2 with unknown mean and variance was recently studied in [START_REF] Di | Statistics for Gaussian Random Fields with Unknown Location and Scale using Lipschitz-Killing Curvatures[END_REF]. The previously cited statistical results permit to derive inference procedures (see, for instance, [START_REF] Lindgren | Spectral moment estimation by means of level crossings[END_REF][START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF][START_REF] Di | Statistics for Gaussian Random Fields with Unknown Location and Scale using Lipschitz-Killing Curvatures[END_REF]) and to test for isotropy, Gaussianity, and marginal symmetry of the underlying fields (see, for instance, [START_REF] Cabaña | Affine Processes: A Test of Isotropy Based on Level Sets[END_REF][START_REF] Berzin | Estimation of local anisotropy based on level sets[END_REF][START_REF] Di Bernardino | A test of Gaussianity based on the Euler characteristic of excursion sets[END_REF][START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF][START_REF] Abaach | Testing marginal symmetry of digital noise images through the perimeter of excursion sets[END_REF][START_REF] Fournier | Identification and isotropy characterization of deformed random fields through excursion sets[END_REF]). [START_REF] Di | Statistics for Gaussian Random Fields with Unknown Location and Scale using Lipschitz-Killing Curvatures[END_REF] propose a test to determine if two images of excursion sets can be compared based on consistent estimators of LKCs. Notice that previous inference methods and tests based on LKCs only require observation of one or two excursion sets, whereas efficient inference about the covariance function and the marginal distribution of X requires more complete observation of the entire field; see [START_REF] Pantle | On the estimation of integrated covariance functions of stationary random fields[END_REF][START_REF] Nieto-Reyes | A random-projection based test of Gaussianity for stationary processes[END_REF], for instance.

Our work provides several important contributions. First, we derive closed formulas for the expected LKCs of Pareto processes with Gaussian or log-Gaussian spectral processes (see Theorem 1). Second, we also show that expected LKCs of mixture processes in the domain of attraction converge to the corresponding expressions of their Pareto process limits when the threshold u used to define excursions in (1) tends to infinity (see Theorem 2). This result is not trivial since no general results exist on the continuity of the excursions LKCs when we consider them as functionals of the stochastic process X. Finally, we derive exact formulas for several specific model classes in non-asymptotic settings.

The paper is organized as follows. Section 2 provides background on the stochastic geometry of excursion sets and presents important assumptions. Formulas for Pareto processes are derived in Section 3 where we directly work with the limit processes without considering convergence towards such processes. The domain-of-attraction setting is treated in Section 4 with some closed analytical formulas of expected LKCs of specific Gaussian mixture constructions and new asymptotic formulas for regularly varying Gaussian mixture processes. Section 5 numerically illustrates convergence rates towards asymptotic expressions. A consistent inference approach for empirical estimation of LKCs using data over regular grids for N = 2 is described in Section 6. We conclude with a discussion and an outlook to follow-up work. Technical lemmas, certain proofs and auxiliary results are postponed to Appendix A. Supplementary materials are provided in Appendix B.

2 Fundamental definitions and assumptions

General Gaussian setting

Throughout the paper, we write W for a Gaussian random field that satisfies the following assumption.

Assumption A0: The process W = {W (s) : s ∈ R N } is a stationary isotropic Gaussian random field on R N , whose sample paths are almost surely in C 3 (R N ). The gradient vector of W at the origin 0 ∈ R N and the Hessian matrix of W at the origin are both non-degenerate multivariate Gaussian variables. Finally, the variable W (0) has zero mean and unit variance, i.e., W is standard Gaussian. Given a nonempty compact domain S ⊂ R N , Assumption A0 ensures that excursion sets in (1) associated to W have positive reach, i.e., that any point outside of S but closer than some fixed positive distance has a unique nearest point in S. We refer the reader interested in more technical background to [START_REF] Thäle | 50 years sets with positive reach -a survey[END_REF] and [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF], for instance.

Remark 1. The stationarity, isotropy and differentiability conditions in A0 imply that

Cov(∂ i W (0), ∂ j W (0)) = λ 2 δ ij , 1 ≤ i, j ≤ N, (3)
where δ ij stands for the Kronecker symbol, and λ 2 is a non-negative constant that is non zero under A0 due to non-degeneracy condition. In the Gaussian literature, the constant λ 2 is called the second spectral moment of W .

In Appendix B, we discuss two commonly used parametric correlation functions whose associated Gaussian field satisfies the above Assumption A0. The correlation function fully characterizes the distribution of a standard Gaussian process, and in the case of a stationary isotropic field it is such that Cor(W (h), W (0)) = σ( h ) with σ(0) = 1, i.e. Cor(W (s + h), W (s)) only depends on the spatial distance h , for any s, h in R N .

Lipschitz-Killing curvatures

We denote by S a fixed compact domain in R N that we choose either as a cube [0, a] N or as a Euclidean ball B N (0, a) with a > 0, if not stated otherwise. Geometric features of exceedance regions are summarized in the Lipschitz-Killing curvatures (LKCs), also called intrinsic volumes or Minkowski functionals in the literature. For any convex compact set D in R N , the LKCs of D, denoted by L j (D) for j = 0, 1, . . . , N , are defined through Steiner's formula (see for instance [START_REF] Robert | Random fields and geometry[END_REF], Section 6.3) that characterizes the volume of the tube D ⊕ B N (0, ρ) with radius ρ ≥ 0, where and| • | is the Lebesgue measure in the corresponding Euclidean space. Then, |D ⊕ B N (0, ρ)| can be expanded as a polynomial function of the radius ρ, whose coefficients are given as follows:

D ⊕ B N (0, ρ) = {x + y : x ∈ D, y ∈ B N (0, ρ)} = {x ∈ R N : dist(x, D) ≤ ρ},
|D ⊕ B N (0, ρ)| = N j=0 ω N -j L j (D) ρ N -j , (4)
with ω j the volume of the j-dimensional unit Euclidean ball,

ω j = |B j (0, 1)| = π j/2 Γ(1 + j/2)
.

For more general compact but not necessarily convex sets in R N , such as the exceedance regions A(u, X, S), one can refer to [START_REF] Robert | Random fields and geometry[END_REF] to get a precise definition of the LKCs. The functional L N is always the N -dimensional Lebesgue measure, as can be seen by taking ρ = 0 in (4). The functional L N -1 is half the (N -1)-dimensional Lebesgue measure of the boundary, and L 0 is the Euler characteristic of D. For the specific domains S ⊂ R N given as a hypercube or a hyperball in R N , we have

L j [0, a] N = N j a j , L j (B N (0, a)) = N j ω N ω N -j a j , j = 0, 1, . . . , N.
For convenience, we write these values more explicitly in Table 3 in Appendix B for the cases of N = 2 and N = 3. For most of the results presented in the following, we could allow for a more general domain S given as any compact stratified submanifold of R N ; see the precise definition in Chapter 8 of [START_REF] Robert | Random fields and geometry[END_REF].

In the remainder of the paper, we use the notation Φ for the standard normal cumulative distribution function, and Φ = 1 -Φ for its survival function.

To conclude this section, we gather the main facts that we need concerning the LKCs of Gaussian exceedance regions in the following proposition. Its first statement is proved, for instance, in [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] Section 3.1. The second result is known as Gaussian Kinematic Formula and can be found in [START_REF] Robert | Random fields and geometry[END_REF] Theorem 13.2.1 or in [START_REF] Adler | Topological complexity of smooth random functions[END_REF] Theorem 4.8.1. Denoting by W (0) the gradient of W at the origin and by W (0) the second derivative (i.e., the Hessian matrix) of W at the origin, one can note that (W (0), W (0)) is non degenerate as required in the quoted theorems, since its two components are independent Gaussian variables which are non degenerate from Assumption A0.

Proposition 1 ([43], [START_REF] Robert | Random fields and geometry[END_REF]). Let W be a Gaussian field on R N satisfying Assumption A0 with positive second spectral moment λ 2 in (3). Then, the following statements hold.

1. There exists a positive constant c such that for any u ∈ R and any k = 0, 1, . . . , N,

E[L k (A(u, W, S)) 2 ] ≤ c |S| 2(1+1/N ) . (5)
2. For any u ∈ R and any k = 0, 1, . . . , N,

E[L k (A(u, W, S))] = N -k j=0 j + k j λ j/2 2 L j+k (S) ρ j (u), (6)
where

j + k j = j + k j ω j+k ω j ω k , and
ρ 0 (u) = Φ(u) , ρ j+1 (u) = -(2π) -1/2 ρ j (u) for j ≥ 0. ( 7
)
Remark 2. The functions (ρ j+1 ) j≥0 defined in ( 7) can be written as

ρ j+1 (u) = (-1) j (2π) -(1+j/2) d j du j (e -u 2 /2 ) = (2π) -(1+j/2) e -u 2 /2 H j (u),
where (H j ) j≥0 are the Hermite polynomials. For later reference, especially for random fields on the geographic space with N = 2, we give the explicit expressions

ρ 1 (u) = (2π) -1 e -u 2 /2 and ρ 2 (u) = (2π) -3/2 ue -u 2 /2 .
Remark 3. The inequality (5) guarantees finite second moments of LKCs. This is an important result for statistical applications since it ensures consistent estimation and asymptotic normality for means of samples from the random variable L k (A(u, W, S)), thanks to the law of large numbers and the central limit theorem, respectively (see Section 6). In specific cases, more precise statements than inequality (5) can be given for the existence of moments of LKCs of excursion sets of a random field. One can heuristically state that a larger index k in {0, 1, . . . , N } corresponds to the existence of higher moments for L k (A(u, W, S)). For instance, for any random field ξ defined on R N , the variable

L N (A(u, ξ, S)) is the N -dimensional volume of the excursion set, such that L N (A(u, ξ, S)) ≤ |S| almost surely. For j = N -1, we have that L N -1 (A(u, ξ, S)) is half the (N -1)-dimensional volume of the level set {x ∈ S : ξ(s) = u}, and it is shown in [7] that E[L N -1 (A(u, W, S)) m ] < ∞ for any m ≥ 0, provided that W is a Gaussian field with almost surely C ∞ sample paths.
On the other hand, the assumption of C 3 sample paths for a Gaussian field W implies the finiteness of the second moment of L 0 (A(u, W, S)), but no general statement can be made about the finiteness of higher moments. The existence of moments of L k (A(u, W, S)) is fundamentally linked with the regularity of W . To the best of our knowledge, so far no general results are available for a precise assessment of this relationship.

Expected LKCs for Pareto processes

In this section, we use the construction of the fundamental point-process limits in functional extremevalue theory to define the class of Pareto processes that arise as certain peaks-over-threshold limits of stochastic processes. We derive formulas for their expected LKCs in the general case, and for the important special cases where Gaussian processes are used in the spectral construction (2) of the corresponding point-process limit. The behaviour of expected LKCs in the domain-of-attraction setting, where we consider the convergence of Gaussian mixtures towards Pareto limit process, will be presented in more detail in Section 4.

Point processes, max-stable processes, Pareto processes

The spectral representation (2) defines a class of Poisson point processes and is the fundamental tool for constructing the different classes of limit processes that can arise for i.i.d. copies of linearly rescaled stochastic processes X. [START_REF] De | A spectral representation for max-stable processes[END_REF] showed that the spectral construction gives a complete characterization of all possible limits when the marginal distribution is prescribed, and he used it to construct the class of simple max-stable processes Z = {Z(s), s ∈ S} (i.e., with unit Fréchet marginal distributions) having continuous sample paths in C(S). Any such max-stable process can be represented as the componentwise maximum of the Poisson process in (2) as follows,

{Z(s), s ∈ S} = ∞ max i=1 Y i (s), s ∈ S . (8)
Recall that {Y i } is a Poisson process whose points Y i are functions on S, and we write M for the intensity measure of this process on C + (S) \ {0}, the set of nonnegative functions excluding the null function. The measure M satisfies the homogeneity property t × M (tB) = M (B), for Borel sets B and t > 0. The tail of the marginal intensity measure of the Poisson process

{Y i (s)} (with s fixed) is M s [z, ∞) = 1/z, z > 0.
Moreover, we can state that sample paths of the points

Y i are in C k if the random function V is in C k .
The full class of possible limit processes is obtained by allowing for marginal transformations

T µ(s),σ(s),ξ(s) (Y i (s)) = µ(s) + σ(s)(Y i (s) ξ(s) -1)/ξ(s)
with deterministic functions for the parameters of location µ(s), scale σ(s) > 0 and shape ξ(s), where the case ξ(s) = 0 is defined as the limit µ(s) + σ(s) log Y i (s). Since the transformation T µ(s),σ(s),ξ(s) is strictly monotonic, the simple max-stable process Z(s) is transformed accordingly to obtain the full class of max-stable processes with margins T µ(s),σ(s),ξ(s) (Z(s)).

The extremal coefficient of S, an important summary of the extremal dependence strength in S, is given by the finite positive value

θ S (M ) = M (B max ) ∈ [1, ∞), where B max = {f ∈ C + (S) : max s∈S f (s) ≥ 1}.
We will often simply write θ S if the structure of M is known from the context. A value of 1 arises for θ S with full spatial dependence in the V i -processes, and larger values correspond to weaker dependence. If interest is not in maxima but in threshold exceedances, then generalized Pareto processes are the only possible limits when conditioning a process X on the exceedance of the process (for at least one location) above a high threshold function ( [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF]). A Pareto process Y always has marginal Pareto tails with shape parameter 1 and scale parameter 1/θ S , and its generalized version allows for more general marginal distributions by applying the same marginal transformations T µ(s),σ(s),ξ(s) (Y ) as before. The Pareto process Y = {Y (s), s ∈ S} associated with the exponent measure M is the stochastic process with probability distribution

P(Y ∈ B) = M (B ∩ B max ) θ S =: M P (B) (9)
for Borel sets B ⊂ C + (S) [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF]. The scaling of the intensity measure of B max by the extremal coefficient θ S in ( 9) ensures that we obtain a probability measure. Pareto processes exhibit Peaks-Over-Threshold (POT) stability, which is characterized through the following equality in distribution, valid for any threshold u ≥ 1:

u -1 Y (•) | max s∈S Y (s) > u d = Y (•).
By construction of the spectral representation (8), any Poisson point Y i falling into B max has distribution according to the Pareto process:

P(Y i ∈ B | Y i ∈ B max ) = M P (B).
We can therefore use the point-process construction (2) to derive properties of the corresponding Pareto process.

Expected LKCs of Pareto processes

The functional point-process representation (2) offers a convenient way to calculate expected Lipschitz-Killing curvatures of the exceedance regions of the Pareto process Y , provided that such expectations are known for the random function V . We here assume that V has sample paths in C 3 (R N ) and is non negative. A threshold u ≥ 1 for Y i (s) corresponds to a threshold u/Λ i for V i (s). Therefore, we can use a conditional expectation argument to obtain the expected LKC of Y by integrating with respect to the intensity function y -2 dy of the Poisson points

{Λ i } in (0, ∞), i.e., if E [L k (A (•, V, S))] is integrable on R + , E [L k (A(u, Y , S))] = 1 θ S ∞ 0 E [L k (A (u/y, V, S))] y -2 dy = 1 u θ S ∞ 0 E [L k (A (y, V, S))] dy, (10)
for k = 0, 1, . . . , N . In the above formula, the leading factor 1/θ S stems from the rescaling of the intensity measure in [START_REF] Berzin | Estimation of local anisotropy based on level sets[END_REF]. The second equation follows from a change of variables from u/y to y. We emphasize that the use of different thresholds corresponds to simple rescalings of the expected LKC. This reflects the peaks-over-threshold stability of Pareto processes. Therefore, we can conveniently define a new quantity that does not depend on u:

C k (Y , S) := E [L k (A(1, Y , S)] = u E [L k (A(u, Y , S)] , u ≥ 1.
We call C k (Y , S) the expected conditional LKC of order k since it can be viewed as a LKC conditional to an exceedance of any threshold u ≥ 1.

Formulas for Gaussian-based Pareto processes

Several specific models given by ( 2) with processes V defined as the transformation of a standard Gaussian process W on C(S) have been studied and used extensively in the extreme-value literature, as already pointed out in the introduction. Recall that a + = max(a, 0) for a ∈ R. A first possibility is to use i.i.d. copies V i of

V d = W α + E[W α + ] , α > 0, (11)
with

E[W α + ] = 2 α/2-1 π -1/2 Γ((α + 1)/2) (12)
in ( 2), leading to the max-stable process Z(s) known as extremal-t process [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF]. The corresponding Pareto processes have been labelled elliptical Pareto processes by [START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF]. A second possibility is to set

V d = e αW E[e αW ] , α > 0, ( 13
)
where E[e αW ] = e α 2 /2 , leading to a large class of max-stable processes of Brown-Resnick type [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]. In both cases, we have E[V ] = 1. We note that the classical Brown-Resnick process as defined by [START_REF] Bruce | Extreme values of independent stochastic processes[END_REF] cannot be handled within our framework since it uses nonstationary processes log V i given by fractional Brownian motion with nonconstant mean function. In the following, we use notation with superscripts scal and loc to refer to the constructions based on ( 11) and ( 13), respectively. The reason for these superscripts will become clear in the next section.

For elliptical Pareto processes, the extremal coefficient has expression [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF]

θ scal S = E [max s∈S W + (s) α ] E W α + . (14)
For Brown-Resnick Pareto processes, the corresponding expression is

θ loc S = E max s∈S e αW (s) E [e αW ] . (15)
Our first important new result states the formulas for expected conditional LKCs of these Gaussian-based Pareto processes.

Theorem 1 (Expected conditional LKCs of Gaussian-based Pareto processes). Let W be a Gaussian random field defined on R N that satisfies Assumption A0.

For the Gaussian-based Pareto processes based on spectral functions in [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF] or (13), the expected conditional LKCs are given as

C k ((Y ) type , S) = c type k (α, λ 2 , S) θ type S , k = 0, 1, . . . , N, type ∈ {scal, loc} (16)
with the corresponding extremal coefficients in ( 14) and (15), and

c type k (α, λ 2 , S) = N -k j=0 j + k j λ j/2 2 L j+k (S) K type j (α), (17)
where the constants K type j (α) are given as follows. For the elliptical Pareto process, we have

K scal 0 (α) = 1, K scal j+1 (α) = α π 1/2 Γ((α + 1)/2) 0≤i≤j β i j 2 i/2 Γ((α + i)/2), (18)
with (β i j ) 0≤i≤j such that ρ j+1 (u) = e -u 2 /2 0≤i≤j β i j u i , for j ≥ 0. For the Pareto process of Brown-Resnick type, we have for j ≥ 0,

K loc j (α) = (2π) -j/2 α j . ( 19
)
Proof of Theorem 1. We obtain the formulas by simplifying the general equality [START_REF] Biermé | On the perimeter of excursion sets of shot noise random fields[END_REF] after inserting the specific Gaussian-based forms of the spectral function V . The detailed proof is postponed to Appendix A.

Remark 4. Unsurprisingly, from Theorem 1 we recover that C N (Y , S) = |S|/θ S for both types of Gaussian-based Pareto processes that we have considered.

Remark 5. Writing out ρ 1 and ρ 2 in explicit form allows one to get

K scal 1 (α) = π -1/2 Γ(1 + α/2) Γ((1 + α)/2) , K scal 2 (α) = α 2π . ( 20
)
The closed-form expressions in [START_REF] Cotsakis | Surface area and volume of excursion sets observed on point cloud based polytopic tessellations[END_REF] for K loc j (α) and in [START_REF] Cotsakis | Statistical properties of a perimeter estimator for spatial excursions observed over regular grids[END_REF] for K scal j (α), j = 1, 2, yield the next results for geographic space with dimension N = 2:

c scal 2 (α, λ 2 , S) = |S|, c scal 1 (α, λ 2 , S) = 1 2 |∂S| + 2 3/2 π -1/2 Γ(1 + α/2) Γ((1 + α)/2) λ 2 |S|, c scal 0 (α, λ 2 , S) = 1 + 1 2 √ π Γ(1 + α/2) Γ((1 + α)/2) λ 2 |∂S| + α 2π λ 2 |S|, c loc 2 (α, λ 2 , S) = |S|, c loc 1 (α, λ 2 , S) = 1 2 |∂S| + 2απ -1/2 λ 2 |S|, c loc 0 (α, λ 2 , S) = 1 + α 2 √ 2π λ 2 |∂S| + α 2 2π λ 2 |S|.

Subasymptotic LKCs for Gaussian mixtures

We now turn to the subasymptotic setting with a focus on Gaussian location or scale mixture processes, which are often used for modeling and simulating in a subaysmptotic framework when the approximation of extreme-value data by asymptotic models is too coarse [START_REF] Huser | Bridging asymptotic independence and dependence in spatial extremes using Gaussian scale mixtures[END_REF][START_REF] Huser | Advances in statistical modeling of spatial extremes[END_REF]. As before, we set a fixed compact domain S in R N , and assume that it is given as a hypercube or hyperball, if not specified otherwise. We will study Gaussian mixture processes {X(s) : s ∈ S} that have the following general structure:

X(s) = g(W (s), Λ), s ∈ S, (21)
where Λ is a shape variable with suitable properties, and with a certain link function g such that the function g(•, λ) is strictly increasing for any λ. This setting allows us to introduce the inverse function h(•, λ) of g so that we have the following almost sure identity for the exceedance region at level u:

A(u, X, S) = A(h(u, Λ), W, S).

Specifically, location mixtures arise for g(w, λ) = w + λ, whereas scale mixtures are obtained by setting g(w, λ) = λw. Under univariate regular-variation conditions expressed in terms of the location or scale variable Λ, the corresponding Gaussian mixture processes converge to the Gaussian-based limit processes (modulo marginal probability integral transformations) previously discussed in Section 3.

In this section, we first derive the expressions of expected LKCs for general Gaussian location or scale mixtures. Then, we use these results to characterize the tail behaviour of expected LKCs in the domainof-attraction setting based on regular variation, and we show that the obtained expressions converge to those derived for the limit processes in Theorem 1. This is an important result in itself since there are no general results on the continuity of LKCs of excursion sets when they are viewed as functionals of the underlying process, such that we cannot invoke any continuous-mapping argument.

Expected LKCs for Gaussian mixture processes

In order to extend results from Gaussian fields (see Section 2) to general Gaussian mixtures X = g(W, Λ) as in [START_REF] Anthony | Statistics of extremes[END_REF], the key point is the conditional expectation representation with respect to the shape variable Λ. It allows us to write, for any function ϕ such that E|ϕ(A(u, X, S))| < ∞,

E[ϕ(A(u, X, S))] = E [E [ϕ (A(h(u, Λ), W, S)) | Λ]] . (22)
We consider a Gaussian field W satisfying Assumption A0 and a random variable Λ independent of W . Let X be the Gaussian mixture X = g(W, Λ) introduced in [START_REF] Anthony | Statistics of extremes[END_REF]. Then, by applying [START_REF] Anthony C Davison | Geostatistics of dependent and asymptotically independent extremes[END_REF] and Proposition 1, there exists a constant c such that, for all u ∈ R and k = 0, 1, . . . , N ,

E[L k (A(u, X, S)) 2 ] ≤ c |S| 2(1+1/N ) ,
and

E[L k (A(u, X, S))] = N -k j=0 j + k j λ j/2 2 L j+k (S) E[ρ j (h(u, Λ))]. ( 23
)
We now focus on expected LKCs for Gaussian location or scale mixtures. We define a Gaussian scale mixture random field by prescribing

X(s) = Λ W (s), s ∈ R N , {W (s)} ⊥ Λ, (24)
where Λ > 0 can be viewed as a random standard deviation parameter embedded in the Gaussian random field W . Early multivariate distributional characterizations (i.e. for a finite number of locations s) were given by [START_REF] Kelker | Infinite divisibility and variance mixtures of the normal distribution[END_REF][START_REF] Cambanis | On the theory of elliptically contoured distributions[END_REF], and an overview is given in [START_REF] Arjun K Gupta | Elliptically contoured models in statistics and portfolio theory[END_REF] for multivariate properties and in [START_REF] Ma | Construction of non-gaussian random fields with any given correlation structure[END_REF][START_REF] Ma | Elliptically contoured random fields in space and time[END_REF] for spatial processes. Constructions with closed-form expressions of multivariate probability densities include the Student-t fields [START_REF] Røislien | T-distributed random fields: a parametric model for heavy-tailed well-log data[END_REF][START_REF] Ma | Student's t vector random fields with power-law and log-law decaying direct and cross covariances[END_REF], Laplace fields [START_REF] Opitz | Modeling asymptotically independent spatial extremes based on Laplace random fields[END_REF] or slash fields [START_REF] Héctor W Gómez | A new family of slash-distributions with elliptical contours[END_REF][START_REF] Huser | Bridging asymptotic independence and dependence in spatial extremes using Gaussian scale mixtures[END_REF], for which the variable Λ 2 follows an inverse gamma, exponential or Pareto distribution, respectively. Similarly, we define the Gaussian location mixtures:

X(s) = Λ + W (s), s ∈ R N , {W (s)} ⊥ Λ, (25)
with Λ a random shape variable. Such processes have received considerable attention in the recent spatial statistics literature due to their ability to account for asymmetric lower and upper tails, and for extremal dependence that is stronger than in the purely Gaussian case (see e.g., [START_REF] Krupskii | Factor copula models for replicated spatial data[END_REF]). In spatial extreme-value theory, they lead to extensively studied limit processes of Brown-Resnick type (see e.g., [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF][START_REF] Jennifer | Efficient inference for spatial extreme value processes associated to log-Gaussian random functions[END_REF][START_REF] Dombry | Exact simulation of max-stable processes[END_REF]).

For a Gaussian scale mixture process, Equation ( 23) applies with function h(u, λ) = u/λ. It is then possible to get closed analytical formulas for the functions E[ρ j (•/Λ)], j = 0, 1, 2, (and associated E[L k (A(u, X, S))] via Equation ( 23)) in the case of specific scale variables (see Table 1). In Table 1, γ(a, •) stands for the lower incomplete Gamma function, i.e., γ(a, x) = x 0 t a-1 e -t dt, for a > 0, x > 0, Γ(a, •) stands for the upper-incomplete Gamma function, i.e., Γ(a, x) = ∞ x t a-1 e -t dt, for a > 0, x > 0, and K k is the modified Bessel function of second kind of order k. Note that Γ

(k, x) = (k -1)! e -x 0≤j≤k-1 1 j! x j , x > 0, for non-negative integers k. Hence, if Λ 2 ∼ Gamma(k, θ), with k ∈ Z, k ≥ 1, then E[ρ 0 (u/Λ)] = 1 2 e -u √ 2/θ u 2/θ 0≤j≤k-1 u 2j j!θ j M u √ 2/θ (1 -j),
where M µ (j) stands for the j-th moment of the Inverse Gaussian distribution IG(µ, µ 2 ).

For a Gaussian location mixture process in [START_REF] De | A spectral representation for max-stable processes[END_REF], Equation ( 23) applies with h(u, λ) = u-λ. By definition of ρ j (see Equation ( 7)), we get the recursion

E[ρ 0 (u -Λ)] = E Φ(u -Λ) , E [ρ j+1 (u -Λ)] = (2π) -1/2 d du E[ρ j (u -Λ)], Law of Λ 2 E[ρ i (u/Λ)] Gamma(k, θ) 1 √ 2πΓ(k) R+ e -(z 2 /2) Γ k, u 2 θz 2 dz, i = 0 with k > 0, θ > 0 2 1-k/2 u k Γ(k)θ k/2 K k ( √ 2θu), i = 1 -(2π) -1 u d du E[ρ 0 (u/Λ)], i = 2 Exp(θ) 1 2 e - √ 2/θu , i = 0 with θ > 0 2/θu K 1 ( 2/θu), i = 1 1 4π 2/θ u e - √ 2/θu , i = 2 Pa(α) u -2α 2 α-1 (π) -1/2 γ(α + 1/2, u 2 /2) + Φ(u), i = 0 with α > 0 u -2α 2 α-1 (π) -1 α γ(α, u 2 /2), i = 1 u -2α 2 α-1 (π) -3/2 α γ(α + 1/2, u 2 /2), i = 2 Table 1: Computations of E[ρ j (•/Λ)], j = 0, 1, 2, for specific distributions of the scale variable Λ 2 . Law of Λ E[ρ 0 (u -Λ)] N (m, σ 2 ) Φ u-m √ 1+σ 2 Gamma(k, θ) 1 √ 2πΓ(k) R e -(z 2 /2) Γ(k, (u -z) + /θ) dz k > 0 and θ > 0 Exp(θ) e 1/(2θ 2 ) e -u/θ Φ(u -1/θ) + Φ(u) θ > 0 Table 2: Computations of E[ρ 0 (u -Λ)]
for specific distributions of the location variable Λ.

for j ≥ 0. Hence, closed-form expressions of E[ρ j (u -Λ)], j ≥ 1, are available provided that we can compute E[ρ 0 (u -Λ)]. In Table 2 we perform such computations for specific distributions of the location variable Λ.

We continue this section with a focus on extreme values of scale and location mixtures by considering regularly varying (RV) scale random variables Λ and exponentially regularly varying location random variable Λ. We refer to Appendix A for definitions and some well known facts on random variables that are regularly varying, i.e., possessing a distribution with power-law tail.

The following result shows the asymptotic behaviour of the expected LKCs of exceedance regions of X = Λ W when Λ is regularly varying at infinity with index α > 0, i.e., Λ ∈ RV ∞ (α) and of X = Λ + W when e Λ is regularly varying in RV ∞ (α) with α > 0, i.e., Λ is exponential-tailed with positive rate α. We use the notation f (u) ∼ g(u) for lim (i) Denote by X scal the Gaussian scale mixture Λ scal W in [START_REF] De | High-dimensional peaks-over-threshold inference[END_REF] where Λ scal ∈ RV ∞ (α) with α > 0 is a positive random variable, independent of W and with continuous probability density. Then, for k = 0, 1, . . . , N , the function

u → E[L k (A(u, X scal , S))] belongs to RV ∞ (α).
More precisely, for k = 0, 1, . . . , N , the following asymptotic holds:

E[L k (A(u, X scal , S))] ∼ u→∞ u -α L Λ (u) 1 c α c scal k (α, λ 2 , S), (26)
where c α := (E W α + ) -1 and L Λ is the slowly varying function given by L Λ (u) = u α P(Λ scal > u). (ii) Denote by X loc the Gaussian location mixture defined by Λ loc + W in Equation [START_REF] De | A spectral representation for max-stable processes[END_REF] where Λ loc is a random variable, independent of W , with continuous probability density such that e Λ ∈ RV ∞ (α) with α > 0. Then, for k = 0, 1, . . . , N , the function

u → E[L k (A(u, X loc , S))
] is exponentially regularly varying with index α.

More precisely, for k = 0, 1, . . . , N , the following asymptotic holds:

E[L k (A(u, X loc , S))] ∼ u→∞ e -αu L Λ (e u ) e α 2 /2 c loc k (α, λ 2 , S), (27)
where L Λ is the SV function given by L Λ (e u ) = e αu P(Λ > u).

Notice that the asymptotics in ( 26) and ( 27) for k = N can both be derived from the well-known Breiman's lemma (see [START_REF] Breiman | On some limit theorems similar to the arc-sin law[END_REF][START_REF] Daren | Subexponentiality of the product of independent random variables[END_REF] and [START_REF] Anthony | Convolution equivalence and infinite divisibility[END_REF], Lemma 2.1). For the sake of clarity, this classical result is recalled in Appendix (see Lemma 1).

Remark 6 (Sub-Gaussian random fields). Theorem 2 in [START_REF] Adler | Excursion sets of three classes of stable random fields[END_REF] is a particular case of our Proposition 2, obtained by focusing on the Euler-Poincaré characteristic (k = 0) for a Gaussian scale mixture with scale variable Λ following the α-stable distribution with α ∈ (0, 1) (see last row of Table 4 in Appendix B). This type of random field represents a simple variant among the stable random fields, called subGaussian in [START_REF] Adler | Excursion sets of three classes of stable random fields[END_REF]. In this case, the function L Λ tends to a known constant depending on α, which allows us to recover the following estimates by applying Proposition 2 (see (4.1) in [START_REF] Adler | Excursion sets of three classes of stable random fields[END_REF]):

E[L 0 (A(u, X, S))] ∼ u→∞ u -α 1 Γ(1 -α) N j=0 λ j/2 2 L j (S) K scal j (α),
with K scal j (α) as in Equation [START_REF] Daren | Subexponentiality of the product of independent random variables[END_REF].

Proof of Proposition 2. In order to get ( 26) and ( 27), we use [START_REF] Anthony C Davison | Statistical modeling of spatial extremes[END_REF] and plug in the respective asymptotics [START_REF] Fournier | Identification and isotropy characterization of deformed random fields through excursion sets[END_REF] and ( 37) that are provided in Lemma 2 in Appendix A.

Based on results summarized in [START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF] (page 59), in Table 4 in Appendix B we gather a list of examples of the distribution of Λ for which we can provide explicit expressions of the asymptotic behaviour of the SV function L Λ .

Domain-of-attraction setting for Pareto processes and expected-LKC convergence

We use the previous results to characterize the tail behaviour of expected LKCs in the domain-ofattraction setting based on regular variation, and we show that the obtained expressions converge to those derived for the limit processes in Theorem 1. The equality of the expressions of LKCs of limit processes and limits of LKCs is not trivial since there are no general results about the continuity of LKCs, such that continuous-mapping arguments cannot be exploited to prove this equality directly.

A natural approach to study peaks-over-threshold limits of a stochastic process X in a compact domain S is to condition on an exceedance of the spatial maximum as detailed by [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF]. If a limit exists then it is a Pareto process, and such a limit exists if and only if the equivalent limits for point processes and componentwise maxima exist. Generalized Pareto processes arise as limits when more general marginal distributions are allowed, but we can always achieve a standard Pareto process limit through marginal pre-transformations of the original process X. Therefore, to study properties related to the extremal dependence of X, we consider the normalized process

X (•) = 1/ F (X(•)), where F (u) = P(X(s) > u). ( 28
)
It possesses standard Pareto marginal distributions with P(X (s) > v) = 1/v, for v ≥ 1, if the marginal distributions of X are continuous, as assumed in our setting. If the following functional convergence in distribution is satisfied as u tends to infinity,

u -1 X (•) | max s∈S X (s) > u d -→ Y (•),
then the limit process Y (•) is a Pareto process [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF]. Its extremal coefficient arises as follows,

θ S = lim u→∞ θ S (u), θ S (u) = u P(max s∈S X (s) > u). ( 29
)
Elliptical Pareto processes arise as limits for Gaussian scale mixtures X(s) = ΛW (s) with Λ ∈ RV(α), α > 0, whereas Gaussian location mixtures X(s) = Λ + W (s) with Λ exponential-tailed possessing rate α > 0 tend to the Pareto limit processes of Brown-Resnick type.

Theorem 2 (Expected conditional LKCs for Gaussian mixture limits). Let W be a Gaussian field defined on R N that satisfies Assumption A0.

Denote by X scal the Gaussian scale mixture X = Λ scal W in [START_REF] De | High-dimensional peaks-over-threshold inference[END_REF] where Λ scal ∈ RV ∞ (α) with α > 0 is a positive random variable, independent of W and with continuous probability density.

Denote by X loc the Gaussian location mixture defined by X = Λ loc + W in Equation ( 25) where Λ loc is a random variable, independent of W , with continuous probability density, such that e Λ ∈ RV ∞ (α) with α > 0.

Finally, let (X ) type and (Y ) type be the corresponding normalized process and Pareto limit process, respectively. Then, for k = 0, 1, . . . , N and type ∈ {scal, loc},

lim u→∞ E L k (A(u, (X ) type , S)) | max s∈S (X ) type (s) > u = C k (Y ) type , S , (30)
with the conditional LKCs on the right-hand side computed in Theorem 1.

Proof of Theorem 2. Let us first remark that the conditional expectation in the left-hand side of ( 30) can be reduced to the ratio E L k (A(u, (X ) type , S)) /P(max s∈S (X ) type (s) > u) since the excursion A(u, (X ) type , S) is empty as soon as max s∈S (X ) type (s) ≤ u.

To exploit the previously established results, we rewrite

E [L k (A(u, X , S))] = E [L k (A(T (u), X, S))],
where T (u) = F -1 (1/u) tends to infinity as u goes to infinity. We are then in position to use the asymptotics given in of Proposition 2.

In the scale mixture case, Equation ( 26) allows us to get that

E [L k (A(T (u), X, S))] ∼ T (u) -α L Λ (T (u)) 1 c α c scal k (α, λ 2 , S).
Let us focus on the term T (u) -α L Λ (T (u)) 1 cα . On the one hand,

P(X(s) > T (u)) = P(W (s) > T (u)/Λ) = E[ρ 0 (T (u)/Λ)] ∼ T (u) -α L Λ (T (u)) 1 c α ,
by using Lemma 2 in Appendix A to get the asymptotics. On the other hand, P(X(s) > T (u)) = P(X (s) > u) = 1/u. Hence, T (u) -α L Λ (T (u)) 1 cα ∼ 1/u. Thanks to [START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF], the equality in ( 30) is established in the Gaussian scale mixture case. In the location mixture case, we can follow a similar procedure to establish the result starting from [START_REF] Di Bernardino | A test of Gaussianity based on the Euler characteristic of excursion sets[END_REF].

Numerical illustrations

In this section we provide graphical illustrations focusing on two-dimensional processes (N = 2). Important technical details of the underlying estimation algorithms for the random quantities L k (A(u, X, S)), based on observations over a regular square grid on R 2 , are discussed later in Section 6. We consider a scale mixture model X = ΛW with a Bergmann-Fock Gaussian random field W with covariance as in Example 1 in Appendix B with a = (100/m) 2 , in a domain of size m × m with m = 2 10 . Sampling is achieved through Matlab using the approach of circulant embedding matrices. We take two specific distributions for the scale variable Λ. In the first row of Figure 1, we consider Λ 2 to be exponentially distributed (Λ 2 ∼ Exp(1)); in the second row, Λ 2 is Pareto distributed with parameter 4 (Λ 2 ∼ Pa(4)). In the second case we have Λ 2 ∈ RV(4) (see Section 4).

In the first numerical study, gathered in Figure 1, we illustrate asymptotics in Proposition 2 and the closed analytical formulas of expected LKCs of Gaussian scale mixtures provided by Table 1. (center panels) and E[L 2 (A(u, X, S))] (right panels) are displayed in logarithmic scale with blue curves (see Table 1). Empirical means of estimated L k (A(u, X, S)) from 500 i.i.d. samples (red stars) are given for several values of u. In the second row, theoretical values of the asymptotic expected LKCs from Equation [START_REF] De Haan | Extreme value theory: an introduction[END_REF] in Proposition 2 are shown as green curves.

Secondly in this section we aim to illustrate the asymptotic behaviour provided by Theorem 2, specifically Equation [START_REF] Dombry | Exact simulation of max-stable processes[END_REF]. We use the same scale mixture model as before but now with different parameters in the two distributions that we consider for Λ 2 : Pareto distribution with parameters α = 3/2 (row (a) in Figures 2 and3) and α = 1/2 (row (b) in Figures 2 and3).

In these two particular mixture cases the marginal distributions of X belong to the Slash family with density functions f (x) = (1-e -x 2 /2 )( √ 2πx 4 ) -1 , for x = 0, and

f (x) = (3(2-(2+x 2 )e -x 2 /2 ))( √ 2πx 4
) -1 , for x = 0, respectively. These Slash marginal distributions are used here to build the standardized process X (see Equation ( 28)).

Numerical illustrations of Equation ( 30) of Theorem 2 are displayed in Figure 2, where the y-axis is given on logarithmic scale. For a third illustration, we write the extremal coefficient θ scal S as in [START_REF] Bruce | Extreme values of independent stochastic processes[END_REF] as

θ ( * ) S := E [max s∈S W + (s) α ] 2 α/2-1 π -1/2 Γ((α + 1)/2) . (31)
Alternatively, we can estimate it using a large value u by one of the following approximations:

θ ( ) S (u) := u P(max s∈S X (s) > u), (32) θ (k)
S (u) :=

c scal k (α, λ 2 , S) E[L k (A(u, X , S))] P(max s∈S X (s) > u), for k = 0, 1, 2. (33)
The estimate in (32) (resp. in [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF]) is based on Equation (29) (resp. on a combination of Equations ( 16) and ( 17)). In these numerical studies the parameters α and λ 2 in (31)-( 33) are supposed to be exactly known. In Figure 3 (left panels) we provide an illustration of the finite-sample behaviour of the empirically estimated quantities in ( 31)-( 33) based on our simulations using several large values of u. Right panels of Figure 3 show the behaviour of the empirically estimated P(max s∈S X (s) > u) for the same simulations and levels u. The mixture models and defining parameters are the same as in Figure 2.

Obviously, the estimation performance of θ

S for k = 0, 1, 2 in Equation ( 33) is strongly determined by the ability to properly estimate L k from observations over a regular square grid on R 2 , as investigated in Section 6 below. S (u) in [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF], for several large values of u and k = 2 (blue crosses), k = 1 (green crosses) k = 0 (red crosses), evaluated using 5000 i.i.d. simulations. Right panels: empirical estimated P(max s∈S X (s) > u), for the same levels u and the same simulations as in the corresponding left panels.

Empirical estimation of moments of LKCs

Theoretical results derived in previous sections describe expectations of LKCs of exceedance regions of Gaussian mixtures (see Section 4), especially for extreme quantile levels with exceedance probabilities tending to zero. We now discuss how such moments can be estimated numerically from datasets with observations over a regular grid in geographical space R 2 , i.e., N = 2. In practical applications, we often are given a discretized representation of the smooth excursion set from a pixelated image. With Gaussian mixtures X as defined in [START_REF] De | High-dimensional peaks-over-threshold inference[END_REF] and [START_REF] De | A spectral representation for max-stable processes[END_REF], the properties related to the shape variable Λ are statistically identifiable only if we have replicated observations X i , i = 1, . . . , n, with n > 1, of the mixture process. If we had only a single observation of X(•) = g(W (•), Λ), we could at most estimate the realized value of Λ, but not make any useful inferences on the probability distribution of Λ. Therefore, we suppose that the processes {X i } are independent and identically distributed according to the marginal distributions of X. This is exactly the framework of our previous numerical studies (see Figures 1,2 and 3).

The estimation issue of L k : pixel images instead of smooth excursions. Inspired by real-life applications, we assume that processes {X i } have been observed over a square regular grid, i.e., we observe {X i (j 1 , j 2 )}, for (j 1 , j 2 ) ∈ G δ := ((δZ) × (δZ)) ∩ S, where S is the observation window -for example a rectangle -and where δ > 0 is the pixel size. In this pixel setting, for a given threshold u ∈ R, the associated replicated excursion set is given by binary black-and-white images representing spatially dependent Bernoulli variables I(X i (j 1 , j 2 ) ≥ u) for (j 1 , j 2 ) ∈ G δ , where we color in black (resp. white) pixels for values equal to one (resp. zero).

In the literature, various empirical approaches and related asymptotic results have been proposed to deal with estimation of Lipschitz-Killing curvatures of excursion sets observed over G δ . Often, the formulas for estimators are based on the implicit assumption that the process X is only partially observed at the centers of the grid cells [δj 1 , δ(j 1 + 1)) × [δj 2 , δ(j 2 + 1)) in G δ , and it has constant value within each grid cell. In this pixelated setting, with increasingly fine grids as δ → 0, the area statistics, obtained by counting relevant pixels in G δ and multiplying by δ 2 , tends to its true value, but more specific adjustments are required to provide the desirable so-called multigrid convergence property (i.e., strong consistency as the pixel size tends to zero) for statistics involving characteristics such as the perimeter or the Euler-Poincaré characteristic.

For instance, several local counting algorithms for perimeter calculation (LCAPs) are known to converge to a biased estimate of the perimeter. Indeed, LCAPs in the Cartesian square grid G δ converge to the perimeter rescaled by a dimension-dependent constant. In the isotropic framework, this normalization factor is 4/π for N = 2, see [START_REF] Miller | Alternative Tilings for Improved Surface Area Estimates by Local Counting Algorithms[END_REF] where also other types of tessellation in dimensions 2 and 3 are considered and associated biases are evaluated. This pixelization error arising in the perimeter estimation of excursion sets of random fields, and ways to correct it, have recently been investigated in [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF][START_REF] Cotsakis | Statistical properties of a perimeter estimator for spatial excursions observed over regular grids[END_REF][START_REF] Cotsakis | Surface area and volume of excursion sets observed on point cloud based polytopic tessellations[END_REF].

A detailed description of estimation approaches for L k would require a highly technical exposition, for which we refer to papers discussing the multigrid convergence property of LKCs and increasing-domain asymptotics when the observation window S increases towards R 2 (see, e.g., [START_REF] Abaach | Testing marginal symmetry of digital noise images through the perimeter of excursion sets[END_REF][START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF][START_REF] Ebner | Goodness-of-fit tests for complete spatial randomness based on Minkowski functionals of binary images[END_REF][START_REF] Di Bernardino | A test of Gaussianity based on the Euler characteristic of excursion sets[END_REF][START_REF] Cotsakis | Statistical properties of a perimeter estimator for spatial excursions observed over regular grids[END_REF]).

Another issue is that many existing estimators of LKCs may suffer from observation bias due to the intersection of the excursion set in (1) with the observation window S. This additional bias of estimates of L k , for k = 0, 1, can be removed by using an edge correction procedure, see [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF], for instance.

Considered L k estimators. Estimation results shown in this paper (see Figures 1, 2 and3) have been obtained with state-of-the-art LKC estimators. To estimate L k (A(u, X i , S)) from {X i (j 1 , j 2 )} with (j 1 , j 2 ) ∈ G δ , we use algorithms implemented in the Matlab functions bweuler (for k = 0) and bwarea (for k = 2). In contrast to the bweuler and bwarea functions, the existing Matlab function for the perimeter (bwperim) performs poorly because of the pixelisation error. Instead, we implement the procedure recently proposed by [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF] and [START_REF] Abaach | Testing marginal symmetry of digital noise images through the perimeter of excursion sets[END_REF]. Our algorithm is summarized in Algorithm 1 below. We denote by L k (i; u, δ), k = 0, 1, 2, the estimator of the kth LKC obtained by following Algorithm 1 for the ith replicate X i observed over a regular square grid G δ with pixel length δ.

Algorithm 1 Estimation algorithm for L k (i; u, δ), k = 0, 1, 2 (case N = 2) (Step 1) We implement L 0 from bweuler, L 2 from bwarea and L 1 recently proposed by [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF] and [START_REF] Abaach | Testing marginal symmetry of digital noise images through the perimeter of excursion sets[END_REF].

When it is required to specify the connectivity of adjacent pixels, we use the average between the 4-connectivty ("rook" scheme) and the 8-connectivity ("queen" scheme).

(Step 2) Since L 1 is the half perimeter, we divide our L 1 in Step 1 by two.

(Step 3) Since our setting is isotropic, we correct the estimate perimeter in Step 2 by the normalization factor 4/π (see Proposition 5 in [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF], see [START_REF] Miller | Alternative Tilings for Improved Surface Area Estimates by Local Counting Algorithms[END_REF] and [START_REF] Cotsakis | Surface area and volume of excursion sets observed on point cloud based polytopic tessellations[END_REF]).

Notice that L0 and L2, estimated on a square lattice, do not need this bias correction.

(Step 4) Finally, we correct the boundary contribution of ∂S in L 0 (from Step 1) and L 1 (from Step 3) by using the edge correction proposed in Proposition 2.5 in [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF].

Notice that L2 does not need this edge correction.

Law of large numbers for L k . Here we do not address increasing-domain asymptotics, but we describe useful asymptotic mechanisms based on the law of large numbers when the number of available replicates n tends to infinity.

We empirically estimate the mth moment of L k (i; u, δ) by 23) and Remark 3, for sufficient conditions. As estimator of the variance of L k (A(u, X, S)), we define σ2

L m k (n; u, δ)) = 1 n n i=1 L k (i; u, δ) m , m = 1, 2, . . . Consistent estimation is possible if the (m + 1)th moment of L k (i; u, δ) is finite, see Propositions 1, Equation (
k (n; u, δ) = L 2 k (n; u, δ)) -L 1 k (n; u, δ) 2 .
Finally, we compute the empirical variance of the expectation estimator L 1 k (n; u, δ)) by using the formula σ2 k (n; u, δ)/n, and the law of large numbers ensures σ2 k (n; u, δ) → 0 as n → ∞. Moreover, if L k (i; u, δ) is a consistent estimator of the true geometric characteristic L k (A(u, X, S)) as δ → 0 (see the previous discussion about multigrid convergence), then we obtain consistent estimation of the expectation of LKCs: lim

δ→0, n→∞ L 1 k (n; u, δ)) = E [L k (A(u, X, S))] .
Estimation for normalized processes. To estimate LKC properties of the normalized process X (s) = 1/(F (X(s))) in (28), we allow for higher flexibility in univariate marginal distributions by assuming that the observed process X obs can have any continuous and not necessarily stationary marginal distribution F s . Therefore, we consider the probability integral transform

X obs (s) = F -1 s (F (X(s))) ,
where the Gaussian mixture process X(s) has marginal distribution F . Then, X (s) = 1/(F s (X obs (s))).

To estimate properties of X , we must either know the marginal distributions F s for all s, or they must be estimated. In simulation experiments, it is possible to work with exactly known marginal distributions F s , but this is usually not the case with real-world data. For inferences on LKCs for a specific level u of X with marginal exceedance probability 1/u (i.e., u is the marginal return period), we have to estimate only the corresponding quantile T s (u) = F -1 s (1-1/u) of the observed process. Estimators Ts (u) 20

of T s (u) must ensure that differentiability properties of sample paths required for estimation of LKC properties are preserved in the standardized process T -1 s (X obs (s)). The quantile Ts (u) (if unknown) can be estimated using an appropriate quantile estimator [START_REF] Daniel P Mcmillen | Quantile regression for spatial data[END_REF]. If T s (u) is estimated consistently, then the continuous mapping theorem ensures consistent estimation of LKCs of X . In the special case where F s does not depend on s, i.e., where we observe the process with stationary margins, we have the alternative possibility to fix the level T (u) for the observed process and then empirically estimate the exceedance probability pu = P(X obs (s) > T (u)) to obtain the normalized threshold û = 1/(1 -pu ). The law of large numbers ensures consistent estimation of u by û with convergence of û to the true value u as the number of replicates n tends to infinity. Note that estimation of p u is equivalent to estimation of the area expectation E [L N (A(u, X, S))], since both are fully determined by F obs (u). Thanks to the continuous mapping theorem, consistent estimation of LKCs is preserved when replacing u by û.

Conclusion

Our results for Pareto limit processes highlight threshold-invariant behaviour of excursion set characteristics conditional to an exceedance due to the property of POT-stability. This property arises asymptotically and is closely related to regular variation of stochastic processes, i.e., RV-properties that hold uniformly over the compact domain S. Such behaviour arises for Gaussian mixture processes, such as location or scale mixtures, and leads to flexibly parametrized limit models. Our theoretical results on moments of geometric summaries of excursion sets of such processes can be viewed as valuable exploratory and inferential tools for spatial extreme-value analysis when data are available on dense and regular spatial grids. They provide a functional perspective on the extremal behaviour of the process, and can be used to study geometric properties of spatial clusters of extreme values. In particular, the question whether POT-stability is present in a real-life phenomenon, leading to spatial extents of clusters that are invariant to the overall magnitude of the event as measured through the maximum over a domain S, can be addressed in a novel way by using information that is only partially captured by the customary bivariate summaries such as bivariate extremal coefficients.

Whereas expected values of LKCs of excursion sets for Gaussian processes satisfying assumption A0 depend on the spatial dependence structure only through the second spectral moment λ 2 , the situation is different for Gaussian mixture processes. For Gaussian-based Pareto processes, the expected LKCs depend additionally on α, and a wide range of functional forms of LKCs in terms of α and λ 2 arises. Finally, in the subasymptotic setting, an even much larger variety of possible forms of LKCs arises and becomes difficult to describe exhaustively. We intend to exploit results presented in this work to build a toolbox of statistical tests and inference procedures for extremal dependencies in spatial data observed on dense grids. This includes climate data obtained through remote sensing or as climate model outputs, such as reanalysis data obtained through conditioning a physical model on observations available from different sources. Another relevant future extension of our methodology is towards space-time models, for instance by considering a 3-dimensional domain S equal to [0, a] 2 × [0, T ] or B 2 (0, a) × [0, T ]. In this space-time setting we will remove the isotropy assumption included in Assumption A0, which is not realistic in that context.

A Proofs and auxiliary results

Proof of Theorem 1. We start with the assessment concerned with elliptical Pareto processes, i.e., where the type is scal.

We write [START_REF] Biermé | On the perimeter of excursion sets of shot noise random fields[END_REF]

In particular, by applying Fubini Theorem, we get

K scal 0 (α) = α c α ∞ 0 e -u 2 /2 √ 2π u 0 z α-1 dz du = c α ∞ 0 e -u 2 /2 √ 2π u α du = 1.
For j ≥ 0, writing ρ j+1 (u) = e -u 2 /2 0≤i≤j β i j u i , we get

K scal j+1 (

B Supplementary materials

Table 3 recalls the formulas of Lipschitz-Killing curvatures of various standard geometric domains in R 2 and R 3 . Table 4 provides information on the tail index and the slowly varying function of a number of commonly used regularly varying distributions. Finally, we present below two commonly used parametric correlation functions where the associated Gaussian field satisfies Assumption A0 in Section 2.

Proposition 2 .

 2 Let W be a Gaussian field defined on R N that satisfies Assumption A0. Let c type k for type ∈ {scal, loc}, be the constants defined in Equation[START_REF] Cambanis | On the theory of elliptically contoured distributions[END_REF] in Theorem 1.

Figure 1 :

 1 Figure 1: Illustration of expected LKCs. Shape variables are given as Λ 2 ∼ Exp(1) (first row) and Λ 2 ∼ Pa(4) (second row). Exact theoretical values of E[L 0 (A(u, X, S))] (left panels), E[L 1 (A(u, X, S))](center panels) and E[L 2 (A(u, X, S))] (right panels) are displayed in logarithmic scale with blue curves (see Table1). Empirical means of estimated L k (A(u, X, S)) from 500 i.i.d. samples (red stars) are given for several values of u. In the second row, theoretical values of the asymptotic expected LKCs from Equation[START_REF] De Haan | Extreme value theory: an introduction[END_REF] in Proposition 2 are shown as green curves.

12 Λ 2 ∼Figure 2 :

 1222 Figure 2: Averaged values (in logarithmic scale) of estimated L k (A(u, X , S)) based on 5000 i.i.d. sample simulations for k = 2 (blue circles), k = 1 (green circles) and k = 0 (red circles) for several large values of u. We also display the theoretical functions c scal k (α, λ 2 , S) u -1 with c scal k in (17) for k = 2 (blue curves), k = 1 (green curves) and k = 0 (red curves).

Figure 3 :

 3 Figure 3: First row: Λ 2 ∼ Pa(3/2). Second row: Λ 2 ∼ Pa(1/2). Numerical illustrations related to extremal coefficients. Left panels: empirical estimation (in logarithmic scale) of θ ( * ) S in (31) (black horizontal line), θ ( ) S (u) in (32) (magenta crosses); θ (k)

1 ECE

 1 introducing the fact that V = c α W α + with c α = k ((Y ) scal , S) = E L k (A(1, (Y ) scal , S)) L k A (y/c α ) 1/α , W, S dy := c scal k (α, λ 2 , S) θ scal S .Using (6), we get that the above numerator is equal toc scal k (α, λ 2 , S) = y/c α ) 1/α dy = α c α ∞ 0 ρ j (z) z α-1 dz.

α) = α c α 0≤i≤j β i j ∞ 0 z= α c α 0≤i≤j β i j 2 α+i 2 - 1 + k j λ j/ 2 2

 02212 i+α-1 e -z 2 /2 dz Γ((α + i)/2), which gives the desired formula in[START_REF] Daren | Subexponentiality of the product of independent random variables[END_REF] by applying[START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF].Let us turn now to the assessment concerned with Pareto processes of Brown-Resnick type. Similar arguments as before yield toC k ((Y ) loc , S) numerator is equal to c loc k (α, λ 2 , S) = N -k j=0 j L j+k (S) K loc j (α),with D = R + in case (i) and D = R in case (ii). Denoting by f Λ the probability density function of Λ, both cases yieldE[ρ j (h(u, Λ))] = R ρ j (h(u, λ)) f Λ (λ) dλ = D ρ j (z) f Λ ( (u, z)) z) F Λ ( (u, z)) dz,where an integration by part has been performed to get the last line. Let us now focus on case (i) and case (ii) separately.(i) If Λ ∈ RV ∞ (α) and if (u, z) = u/z, then for any fixed z ∈ D = R + , F Λ ( (u, z)) = P(Λ > u/z) ∼ u→∞ z α u -α L Λ (u).Hence, applying Lebesgue's Theorem in the last above integral, we haveE[ρ j (h(u, Λ))] ∼ u→∞ u -α L Λ (u) (-1) R + ρ j (z) z α dz,where (-1)R + ρ j (z) z α dz = R + ρ j (z) αz α-1 dz = (c α ) -1 K scal j (α),thanks to[START_REF] Estrade | A Central Limit Theorem for the Euler characteristic of a Gaussian excursion set[END_REF].(ii) If e Λ ∈ RV ∞ (α) and if (u, z) = u -z, then for any fixed z ∈ R,F Λ ( (u, z)) = P(Λ > u -z) ∼ u→∞ e αz e -αu L Λ (e u ).Applying the same arguments as in case (i), we getE[ρ j (h(u, Λ))] ∼ u→∞ e -αu L Λ (e u ) (-1) R ρ j (z) e αz dz, where (-1) R ρ j (z) e αz dz = R ρ j (z) αe αz dz = e α 2 /2 K loc j (α), thanks to (35).

Acknowledgments

This work has been partially supported by the French government through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR-19-P3IA-0002) and through the project ANR MISTIC (ANR-19-CE40-0005).

with

In the same manner, applying again Fubini Theorem yields K loc 0 (α) = 1 and for j ≥ 0, writing ρ j+1 (u) = -(2π) -1/2 ρ j (u) (see Equation [START_REF] Armentano | Conditions for the finiteness of the moments of the volume of level sets[END_REF]) and performing an integration by parts, we get K loc j+1 (α) = α (2π) -1/2 K loc j (α).

A simple induction yields formula in [START_REF] Cotsakis | Surface area and volume of excursion sets observed on point cloud based polytopic tessellations[END_REF].

Definition 1 (Regular variation RV ∞ (α)). We say that a function G is regularly varying at infinity with index α ≥ 0

For a random variable X, we use the shorthand notation X ∈ RV ∞ (α) if the survival function F of X is regularly varying at infinity with index α, i.e., F ∈ RV ∞ (α).

We have F ∈ RV ∞ (α) if and only if F (x) = L(x)x -α with a slowly varying (SV) function L such that L(tx)/L(t) → 1, t → ∞ (i.e., L ∈ RV 0 ).

We recall below the well-known Breiman's lemma (see [START_REF] Breiman | On some limit theorems similar to the arc-sin law[END_REF][START_REF] Daren | Subexponentiality of the product of independent random variables[END_REF] and [START_REF] Anthony | Convolution equivalence and infinite divisibility[END_REF], Lemma 2.1).

Lemma 1 (Breiman's lemma). Suppose X and Y are independent random variables. If X ∈ RV ∞ (α) with α ≥ 0 and Y ≥ 0 with E(Y α+ε ) < ∞ for some ε > 0, then

Lemma 2. Let Λ be a positive random variable with continuous probability density function.

then the following asymptotics holds for u → ∞:

, then the following asymptotics holds for u → ∞:

Proof of Lemma 2. Let us first remark that the desired asymptotics concern E[ρ j (h(u, Λ))] where h(u, λ) equals u/λ in case (i) and respectively equals u -λ in case (ii). We introduce the inverse (decreasing) function such that for all u > 0 and λ > 0,

values as for instance α = 1/2 0 < α < 1 (Lévy distribution) Table 4: A list of classical RV(α) distributions given by their survival function F Λ with associated SV function L Λ .

Example 1. The standard Bergmann-Fock Gaussian field W is characterized by its correlation function σ( h ; a) = e -a h 2 at spatial distance h with a parameter a > 0 related to spatial range, whose inverse 1/a is called the scale parameter; this model is also called the Gaussian correlation function. For this example, λ 2 = 2 a.

Example 2. As second example, we consider the Matérn correlation function given by σ( h ; ν, a) =

, where K ν is the modified Bessel function of the second kind of order ν, and a > 0 is again a parameter related to spatial range. The smoothness parameter ν > 0 defines the Hausdorff dimension of the graph Gr W = {(s, W (s)) : s ∈ R N } and determines the differentiability of the sample paths. The associated second spectral moment is given by λ 2 = a 2 2(ν-1) , and its existence requires ν > 1. Furthermore, for a positive integer k, the sample paths of a Gaussian Matérn field are k times continuously differentiable if and only if ν > k. Here we require the stronger assumption that W is almost surely C 3 , which implies ν > 3.