{'fraction_non_alphanumeric': 0.03280332681017613, 'fraction_numerical': 0.015019569471624266, 'mean_word_length': 4.411834789515488, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inference for extreme spatial temperature events in a changing climate with application to Ireland

We investigate the changing nature of the frequency, magnitude and spatial extent of extreme temperature in Ireland from 1960 to 2019. We develop an extreme value model that captures spatial and temporal non-stationarity in extreme daily maximum temperature data. We model the tails of the marginal variables using the generalised Pareto distribution and the spatial dependence of extreme events by a semi-parametric Brown-Resnick r-generalised Pareto process, with parameters of each model allowed to change over time. We use weather station observations for modelling extreme events since data from climate models involves abstraction and can over-smooth these events. However, climate models do provide valuable information about the detailed physiography over Ireland. We propose novel methods which exploit the climate model data to overcome issues linked to the sparse and biased sampling of the observations. Our analysis identifies a substantial temporal change in the marginal behaviour, but not the dependence structure, of extreme temperature events over the study domain. We illustrate how these characteristics result in an increased spatial coverage of the events that exceed critical temperatures.


Introduction
The Intergovernmental Panel on Climate Change reports an observable change in extreme weather and climate events since around 1950 (IPCC, 2021). Characterisation of extreme temperature events, in particular, is crucial for societal development, for estimating risks, and to enable the mitigation of their effects for many sectors. Common application areas where the risks appear raised include healthcare, economic growth, agricultural disruption, and infrastructure. Brown et al. (2008) observed a warming of both maximum and minimum temperatures since 1950 for most regions indicating an increasing number of warm days, longer heatwaves, and fewer cold extremes.
In Ireland changing extreme temperature behaviour has also been observed. McElwain and Sweeney (2007) found that a warming of both maximum and minimum temperature observations occurred for all sites over 1961-2005. O'Sullivan et al. (2020 showed that the frequency of extreme temperature events for Co. Dublin has increased over the period 1981-2010. Both these approaches considered only the marginal behaviour of extremes. To the best of our knowledge, the only modelling of spatial extreme temperature events in Ireland is by Huser and Wadsworth (2020). They used a gridded Irish temperature data set, which has the potential to be over-smooth relative to the observed process, and fitted this under the assumption of stationarity over time and space.
Here we are interested in developing a model which captures the temporal evolution of spatial extreme temperature events over Ireland. This involves modelling how the marginal distributions vary over space, accounting for spatial dependence within extreme events, and modelling how these two elements vary over time and established measures of climate change. As we are concerned with events in the tails of the distribution, we focus on the extreme data to avoid biases (Davison et al., 2012).
Observational extreme event data are sparse and so they need to be used efficiently. The traditional statistical approach is to model these data with powerful probabilistic characterisations from extreme value theory. This theory provides a parsimonious asymptotic justification for extrapolation which enables us to describe the properties and behaviour of events which are more extreme than those previously observed. The theory by itself however will not account for spatial mapping over heterogeneous geography when the characteristics of the complex spatial events change over time. Here we propose a novel approach to address these issues which exploits the physical knowledge of the climate processes from information given by fine-scale climate model data. We outline our strategies for each approach below.
The theory of univariate extreme values for stationary processes (Leadbetter et al., 1983) and the associated statistical models (Coles, 2001) fully determine a simple, parametric distributional family, namely the generalised Pareto distribution (GPD), as the non-degenerate limit distribution for the normalised excesses of a threshold, as that threshold tends to the upper endpoint of the marginal distribution. We model the marginal distribution tails using the GPD with its parameters changing smoothly over space and time (Chavez-Demoulin and Davison, 2005, Youngman, 2020. Traditional approaches to spatial extreme modelling use max-stable process models (de Haan and Ferreira, 2006), a natural extension of univariate block maximum limit theory. Brown and Resnick (1977) introduced a widely used subclass of these models, derived from Gaussian random fields, known as Brown-Resnick processes (Kabluchko et al., 2009, Padoan et al., 2010, Engelke et al., 2011. These models cannot capture the spatial patterns of observed extreme events and so inference using them can lead to over-estimation of dependence (Huser and Wadsworth, 2020).
Recent developments for modelling spatial threshold exceedances include generalised r-Pareto processes (Dombry andRibatet, 2015, de Fondeville andDavison, 2018) and spatial conditional extremes models (Wadsworth and Tawn, 2019). The former is restricted to processes which exhibit a strong form of dependence only, known as asymptotic dependence, between all sites in the region of interest, whereas the latter also covers asymptotically independent cases (Coles et al., 1999). In Section 4 we define these properties precisely and provide evidence that generalised r-Pareto processes are suitable for daily maximum temperatures over Ireland. We identify extreme spatial fields, based on observations at d sites, as those which exceed a sufficiently high threshold for a risk function r : R d → R + . We model these fields as realisations of a Brown-Resnick Pareto process, which is closed under marginalisation when r is chosen carefully.
Although we have observational daily maximum data for a network of Irish temperature stations, only 33% of these have more than 30 years of data due to differential operational periods and quality controls. The network exhibits selection bias, with the longest records typically being near the coastline. We also have a rich data set of daily maximum temperatures on a fine-grid over Ireland that were generated from a climate model. These data have no missing values or location biases, and incorporate the impact of known geophysical structures on the temperature process.
Climate models involve some abstraction of the processes they model and tend to under-predict extreme events. Various attempts have been made to downscale the climate data to produce a proxy for the observed data which gives a spatial and temporally complete data set. We prefer to let the observational data stand for themselves, particularly in relation to the information they provide about temporal non-stationarity. Our novel method for inference leverages core information from the climate model data, i.e., parameter estimates (within sample quantiles and GPD parameters) over space.
Our paper is organised as follows. Section 2 details the observational and climate model data used. Sections 3 and 4 describe the marginal and dependence modelling of the process respectively, each setting out how we account for their changing behaviour over time. In Section 5 we use the model to explore how the properties of spatial extreme events have changed over time, which is followed by a discussion in Section 6.

Data
Our observational daily temperature data come from three sources giving a network of 126 Irish temperature stations shown in Figure 1 (left). For the Republic of Ireland, data for 51 stations were from MetÉireann's archive and a further 44 sites were extracted from the database of the Copernicus Climate Change Service (C3S); a consortium who collate in situ observations. Data for 31 Northern Ireland sites were obtained from the Met Office Integrated Data Archive System (MIDAS), through the CEDA archive. Collectively these data exhibit large amounts of missing values, with 61% of daily values observed in the last 30 years, and only 2.4% observed before 1960. No single day has data for every station. The average span of data for each station is about 25 years. The sites with the most data tend to be located near the coast (see Figure 1 (left), reflecting that stations were primarily selected for weather forecasting purposes. Climate models represent the best understanding of the physical processes driving weather and climate. We use a regional climate model's output to provide information on daily temperature data on a fine-grid over Ireland. We use the MOHC-HadREM3-GA7-05 model output, when driven by the global climate model MOHC-HadGEM2-ES. These data reflect observed changes in climate over 1951-2005 on a 0.11 2 degree resolution with 558 grid points shown in Figure 1 (right), with the figure illustrating daily temperatures for a generated event which has an extreme average temperature over Ireland.
To model temporal non-stationarity of extreme temperature data it is common to use time alone as a covariate. Additionally, we include time varying causal covariates, in the form of global mean temperature anomalies from the HadCRUT5 gridded data set and CO 2 emissions for Ireland (from Our World in Data; Ritchie and Roser, 2020). Figure 2 plots these covariates over 1950-2019. We did also consider global CO 2 emissions as a covariate also but found this less effective than the Irish emissions; the key difference between these series being the larger year to year variations and the post 2000 decline in the Irish emissions. We denote the climate change covariates by z t = (GMT t , CO2 t ), and note that they are constant over Ireland in a given year. 3 Marginal models

Overview
Let X o (t, s) denote the observational process of maximum daily temperature at site s ∈ S ⊂ R 2 at time t ∈ T ⊂ Z, and let X c (t, s) be the equivalent process from the climate model. Here S denotes Ireland, s corresponds to the vector of latitude and longitude, and we have data on the two processes at S o ⊂ S and S c ⊂ S. We use the subscripts of the processes to identify their type, though when defining notation for methods which apply to both processes we drop the subscripts.
For modelling X(t, s) we account separately for the behaviour below and above a threshold u(s), which is fixed over time. In Section 3.2 we use spatial quantile regression for the data below the threshold and, Section 3.3, for data above the threshold we use the generalised Pareto distribution (GPD) parametric model. Based on the use of standard GPD threshold selection methods (Coles, 2001), we take u(s) as our estimate of the 95% quantile of X at s. We incorporate changes over time in the marginal distribution function only above the threshold, as values lower than this are less critical here (see Section 5). A novelty in our approach comes from using estimates from the X c process to infer features of the X o process given that |S o | |S c |. Specifically we need to identify thresholds, as well as estimating threshold exceedance probabilities and the GPD parameters, for the tails of both the X c and X o processes.

Modelling the body of the distribution
Following the approach of Yu and Moyeed (2001) and the R package evgam (Youngman, 2020) we use the asymmetric Laplacian distribution (ALD) for quantile regression to estimate a range of spatially varying τ th quantiles, u(s; τ ), up to and including the threshold u(s) = u(s; 0.95), i.e., 0 ≤ τ ≤ 0.95. The density function of the ALD τ is where ρ τ (z) = {τ − I(z < 0)}z is the check function, u ∈ R and ψ > 0. We assume that the parameters u and ψ vary smoothly over the spatial domain S and fit a GAM with thin-plate smoothing splines (Wood, 2006) and an ALD likelihood. In all cases we take the package's default selection method for selecting the number basis functions.
For each τ , we estimate u c (s, τ ) from the climate model data via where u c (s; τ ) and log{ψ c (s; τ )} are 2-D spline functions over S for each fixed τ , i.e., with each function of the form β 0 + K k=1 D k d=1 β kd b kd (s) for basis functions b kd and basis coefficients β kd where K is the number of terms in the additive sum (equal to the number of covariates for non-interacting models) and D k the number of basis functions for additive term k. We repeat this inferential approach to estimate u o (s; τ ), including as covariates both s and u c (s; τ ). This allows us to incorporate richer spatial information that is not captured in the observational data set due to the spatial and temporal sparsity of the observational data. Specifically we take where u o and log ψ o are 2-D spline functions for each fixed τ . To provide estimates for all τ we fit this model over a grid of τ values and then use a cubic spline for each s to give a continuous estimate over 0 ≤ τ ≤ 0.95 for u c (s, τ ) and u o (s, τ ). We keep the grid of τ values relatively coarse to avoid issues of quantile estimates crossing. We found that using u c (s; τ ) as the only covariate for τ < 0.95 and only s for τ = 0.95 worked best.

Modelling the tails of the distribution
We assume that the distribution of exceedances of the threshold u(s) follows a generalised Pareto distribution (GPD), see Pickands (1975) and Coles (2001) for the probabilistic justification and properties of this model. The GPD(σ, ξ) has distribution function for y > 0, which is defined in terms of a shape parameter ξ ∈ R and a scale parameter σ > 0 and where we use the notation x + = max(x, 0). When X(t, s) > u(s) the threshold excess is given by with the shape parameter taken to be constant over time, as is typical in GPD modelling. We model the structure of the GPD parameters differently for the climate model and observational data. These structures are detailed below.
When modelling the climate model data we believe we have sufficient observations to fit the model separately at each site s ∈ S c . We do not attempt to impose any connections between the GPD parameters over s ∈ S as the data are spatially consistent by their generation and so it would be wrong to smooth the parameters spatially if we want to capture as much as possible of the relevant geophysical features of the climate model. To capture the temporal non-stationarity in the data we allow σ c to vary smoothly with time. As |S c | = 558 any attempt to fit this model jointly across sites would be computationally intensive, so we fit this model separately site-by-site. A consequence of this choice is that we do not impose spatial smoothness of either the temporal changes in σ c nor in the values of ξ(s). Thus we fit the model where log σ c (t, s) is a 1-D spline in t for each s ∈ S c .
From this model we take {σ c (t, s}; s ∈ S o } as covariates for inference using the observational tail data, and {σ c (t, s}; s ∈ S c } for interpolating that model over space. As we want to allow the temporal trends in the data to stand for themselves as much as possible we consider these covariates at a fixed time t * , which we take as the year 2000. We do this to ensure that we learn all temporal non-stationarities from the observational data and only geographic information is taken from the climate model. Next, we define and model the extreme data from the observational data set Y o (t, s). As the observational data has, on average, fewer observations per site than the climate model data, and has periods of missing data within longer data spans, we take ξ o (s) = ξ o , for an unknown constant ξ o . This choice reduces the risk of having parameter identifiability problems and was supported by extensive exploratory analysis. This approach is common practice in environmental extreme value modelling.
We investigate the effect of incorporating σ c , from expression (6), and the climate change covariates z t = (GMT t , CO2 t ) into the GPD model of the extreme observational data Y o (t, s). The generic form of each of the models we consider is For log σ o we explore GAM models with linear combinations of spline and linear functions of the covariates and their associated interaction terms. The most advanced model described below (denoted Model 5) is the richest model we attempt to fit; the other models are presented for comparison purposes. For all these key models we give their structure as well as AIC and BIC values: Model 1: Here σ o is a 2-D spline of s and 3-D spline of s and t.
Model 2: As Model 1 but with the inclusion of an additive linear term incorporating the covariate σ c (t * , s) from the climate model output. By fixing a time varying function at a specific value t * (we take t * = 2000) we remove the temporal information from the climate model output term. Model 3: As Model 2 but with the addition of GMT t covariate in a 3-D spline term with s. This allows us to see the relative ability of t and GMT t in explaining the temporal non-stationarity.
Model 4: As Model 3 but an additional 3-D spline term for CO2 t and s. Model 5: As Model 4 but with additional linear terms for each of t and CO2 t and a 1-D spline for GMT t , i.e., log σ o (t, s) = β 0 + β 1 t + β 2 CO2 t + β 3 σ c (s, t * ) + p(GMT t ) + p(t, s) + p(GMT t , s) + p(CO2 t , s) (8) where we use p as a generic spline term. We fit these models using pseudo-likelihood, assuming independence between sites and over time. The AIC and BIC values under this likelihood are given in Table 1 we recognise that formally these should be adjusted for use with the pseudo-likelihood, as suggested by Davison et al. (2012), but as we are only looking for substantial changes in value we do not focus on formal comparisons between them. Model 2 gives a comparable fit to Model 1, but has the benefit of offering spatial interpolation to non-observed sites across Ireland. Model 3 offers a substantial improvement by including GMT t (in addition to t, whilst the additional inclusion of CO2 t emissions for Ireland in Model 4 makes a moderate improvement. Model 5 gives a further very substantial improvement through the inclusion of marginal effects for t, GMT t and CO2 t . Across Models 1-5 the GPD shape parameter estimates are almost unchanged, all being negative giving evidence of a finite upper endpoint to the maximum possible temperature at each site in Ireland. Figure 3 provides the spatial estimates of threshold u o (s) and GPD scale parameter for the year 2019, the latter obtained using Model 5. Each plot is over s ∈ S, not simply for the sites s ∈ S o which were used for this model fit. The threshold (Figure 3, left) shows cooler temperature quantiles on the west of Ireland and coastal regions on the south and north coasts. In contrast the scale parameter (Figure 3, right) shows that the most variable excess distribution is on the west coast, with a decay in values from west to east. For the scale parameter we find a substantial benefit from including information from the climate model data, as the estimate of the coefficient β 3 , of σ c , is close to 1, so we are able to infuse into the spatial interpolation key features from the climate model.  shows the variation of the scale parameter σ o over time (averaged over Ireland as a whole and at sites at the four extremities together with a central site). The temporal signal is similar across all six series, with an additive difference over series. The rapid variations over time may at first seem surprising given our use of splines over t. However, as t changes so do both GMT t and CO2 t ; neither of these covariates is smooth over time, although over the sampling period both increase with t (see Figure 2). The improvement in fit through using these non-smooth features with t shows the complications of modelling non-stationarity as inter-annual variations in climate covariates are important to describe observed temporal variations. Despite this roughness over time, it is interesting that we find an increase of up to 50% in the GPD scale parameter from 1960 to the present day, although prior to 1960 we have what appears to be a decline. We investigated this feature by looking at how many records were available; we have more than four times the data in the 1960s than the 1950s, from 8 and 39 stations in the respective decades. Furthermore, all stations (except one) with data from 1950-59 were coastal, which may also be leading to spatial bias in the data source. We suggest that the decline before 1960 is actually an increase in uncertainty due to the lack of available data before that year, though for brevity we do not plot the error bars in Figure 4. Now we use the selected GPD model for the excesses of u(s) to give a model for the marginal distribution function F t,s of X(t, s). Specifically for y > u(s) we have where λ(t, s) = Pr{X(t, s) > u(s)} is the threshold exceedance probability, and H is given by Equation (5). If there is no temporal non-stationarity in X o (t, s) then by construction of the threshold u c (s) = u c (s; τ ), and we would have λ o (t, s) = 1 − τ = 0.05 across s ∈ S. We use a GAM fit to the number of exceedances per day and per site, i.e., a Binomial model with logit{λ o (t, s)} taking the structural form of the covariate model (8). Figure 4 (right) shows estimates of λ o over time for the same sites and spatial averages as in Figure 4 (left). There is a strong similarity in estimated patterns of temporal variations for λ o and σ o . They show a doubling of exceedance probability over the period 1960-2019, with unreliable estimates in the 1950s, as before.

Spatial models 4.1 Standardising data
When modelling dependence between variables with differing marginal distributions and covariates, it is common to first standardise the marginal variables so that they have an identical distribution over variables and covariates (Coles, 2001, de Haan andFerreira, 2006). Here we transform to (unit) Pareto distributions, X P , X P c and X P o , using the same notation as in Section 3. The choice of Pareto marginal scale is ideal for studying asymptotically dependent variables such as r-Pareto processes (de Haan and Ferreira, 2006), but less ideal for asymptotically independent variables, where shortertailed Exponential or Laplace variables are favoured (Wadsworth and Tawn, 2019). We achieve the transformation through use of the probability integral transform, i.e., , for all s ∈ S and all t where the marginal distribution function F t,s takes a different estimated form below and above u(s) as described in Sections 3.2 and 3.3 respectively. Thus Pr(X P (t, s) > y) = y −1 for all y > 1, t and s. To transform from standard Pareto margins back to the original distribution at s and t, then the inverse of Equation (9) can be used.

Classification of extremal dependence type
We now explore the nature of the extremal spatial dependence structure in the processes X P c and X P o , and so omit the temporal dimension of these spatial processes. Following Coles et al. (1999) we estimate the pairwise coefficient of asymptotic dependence, χ. Specifically, for the process X P at sites s i and s j , χ = χ P (s i , s j ) is defined by where the second equality holds given the marginal distribution of X P . If χ P (s i , s j ) > 0 (or equals 0) then process X P is asymptotically dependent (or asymptotically independent) at this pair of sites respectively. The value of χ (0 < χ ≤ 1) determines the degree of asymptotic dependence, with χ increasing as this dependence strengthens. The selection of which type of extremal dependence to use depends on whether or not the process is asymptotically dependent for all s i , s j ∈ S. The tool that is typically used to identify asymptotic dependence, at a pair of sites, is χ P (·, ·; p), where with v p the pth marginal quantile of X P . The behaviour χ P (·, ·; p) is studied as p increases to 1 to see if it decays to zero or stabilises at a non-zero limit. To estimate χ P (s i , s j ) we take v = v p , the pth quantile of X P and use the empirical value of probability (11), by exploiting the replication of the variables over t. We denote this estimator byχ P (s i , s j ; p). We evaluateχ P (s i , s j ; p) for a range of large p (0.95, 0.975, 0.99). We plotχ P (s i , s j ; p) against the Euclidean distance between sites, i.e., h ij = s i − s j . If the process was both stationary and isotropic and we could perfectly estimate χ P for all pairs in a continuous set S, then this plot would be a smooth curve. Instead as S is discrete and the idealised properties do not hold, we expect variation in estimates for identical values of the distance h and so a cloud of points over h. To summarise this cloud we use conditional quantile regression (Youngman, 2020) for a range of quantiles. Figure 5 shows these estimates for both X P c and X P o processes. Despite the climate model having a much richer set of pairs of sites, associated distances, and longer simultaneous data for inference, both processes provide very similar findings. Naturally χ P (s i , s j ; p) decreases with distance but in all cases it is far from zero, even at the longest distances of pairs of sites from S. For X P o we do see more variation in the conditional quantiles than for X P c , with this being most obvious for the short distances. Furthermore, we have thatχ P o (s i , s j ; p) ≤χ P c (s i , s j ; p) for all distances between pairs, suggesting that the X c data is overestimating the extremal dependence in X o .
Most critical for identifying whether the processes are asymptotically dependent is to see if the estimates change with p. Focusing on the conditional median as p increases, we can see a small decline, at all distances. However, even at p = 0.99 these estimates are far from zero for both X P c and X P o . So it seems reasonable to conclude that maximum daily temperature data are consistent with asymptotic dependence over Ireland.

r-Pareto Processes
We now model the extreme values of the process X P (t, s) over s ∈ S, with Pareto distributed marginal variables. We look at each spatial field separately for each t and so to simplify notation we drop the argument t. First we define what we mean by a spatially extreme event, as unlike in univariate extremes there is no natural ordering of multivariate or spatial processes. Here the level of extremity of the stochastic process X P (t, s) is determined by a risk function r(X P ) := r{X P (s) : s ∈ S} ∈ [0, ∞), where the only other constraint on r is that it is homogeneous of order 1, i.e., for any constant c > 0, then r(cx) = cr(x) for any x with non-negative components. de Fondeville and Davison (2018) suggests that r could be taken as the magnitude at one particular site of interest, or the spatial mean, median, maximum or minimum over S. Dombry and Ribatet (2015) show that under weak conditions on the process X P that Pr v −1 X P (s : s ∈ S) ∈ · | r(X P (s : s ∈ S)) > v → Pr{Y P r (s : s ∈ S) ∈ ·} (12) as v → ∞, where {Y P r (s) : s ∈ S} is marginally non-degenerate in any margins. Here Y P r is called the r-Pareto process. If the limit provided in Equation (12) is a good approximation for large v, then those spatial events with a risk function exceeding v will be well-approximated by an r-Pareto process. Crucially, Y P r can be decomposed into two independent stochastic components as follows: where R is Pareto distributed and can be interpreted as the risk of the process, and {W (s) : s ∈ S} is a stochastic process which describes the spatial profile of the extreme event, i.e., the proportion of the risk at each site. By construction we have that R = r(Y P r (s : s ∈ S)) and r(W (s : s ∈ S)) = 1. This characterisation is powerful as the extrapolation to larger events than those observed is determined by R, which has a known distribution, and that the spatial profiles of larger events have exactly the same stochastic properties as less extreme events. We follow Engelke et al. (2015), and de Fondeville and Davison (2018) by modelling W as a log-Gaussian stochastic process which is determined solely by a variogram. We use the Matérn variogram family for inter-site distance h ≥ 0, K ν a modified Bessel function of the second kind, and the postive parameters (α, φ, ν) determine the variance, range and smoothness respectively. We have not found any spatial extremes applications where missing data has been discussed. When encountering missing data it is tempting to remove all observations at that time point from all sites in the network. However, we observed in Section 2 that this tactic would leave us with no data! Fortunately, thanks to the properties of the log-Gaussian process, it is possible to show that the model is closed under marginalisation (Engelke and Hitz, 2020). With missing data we need to be much more careful in selecting a suitable risk function r. The choice of risk function needs to be invariant to the changing dimension of partially observed events, whatever their missing patterns. Hence in our statistical inference at time t, for all t ∈ T , the set of observation times, we take where I o (t, s i ) is the indicator variable for whether X o (t, s i ) is observed or not, d is the number of sampling locations selected for the risk function evaluation using sites s 1 , . . . , s d . Thus the risk function r t is the average of standardised variables over the d sites which were observed at time t. We pick the set of d sites as the d longest recording stations with the least missing data. We allowed for the variance and range parameters of the Matérn variogram (14) to vary over time t and z t while keeping ν constant. We considered a range of constant and log-linear models but find no statistical significant evidence for change in these parameters with respect to the covariates, although such a change was found before we fully tuned our marginal model. The estimated variogram parameters for our best fitting model, under spatial stationarity and time-invariance, determine an estimate of χ P o ; which is presented in Figure 6. Here the risk threshold we used to model extreme spatial events is the 90% quantile of the risks calculated for all observed events. We explored different threshold choices and selected the highest level we could whilst making the usual bias/variance tradeoff for tail selection. On this figure we also plot the pairwise empirical estimator (11) with p = 0.95, this corresponds to a similar plot in Figure 5 (bottom left). The agreement of the two estimators is good, with the fitted model central to the empirical estimates. We do not believe all the variation in the points is noise from estimation, and instead that different distance metric, which accounts for non-stationarity, may be more effective. See Section 6 for discussion of this. Nevertheless, the results we present in Section 5 are for the spatially stationary r-Pareto process.

Temporal changes in spatial extreme events
We present a range of summaries of how spatial daily maximum temperature extreme events in Ireland are changing over the period 1960-2019. First we look at the changes in the marginal quantiles. Figure 7 plots estimates of the level exceeded by daily values with probability 1/36500, i.e., what would be a 100 year level if the process was stationary in time. For simplicity we refer to these as the 100 year levels which change over time. The estimates are based on Model 5. We show estimates for 2019 (left) and the difference between the return level estimates for 2019 and 1960 (right). In the latter a positive value represents a warming of temperatures. The 100 year return level has increased between 1 − 3 • C, with the smallest change in the north-west and the largest over central Ireland. As the GPD shape parameter is quite negative (−0.25), the finite upper endpoints of the marginal distribution change in a similar way. These changes in extreme temperatures over the observed record are substantially larger than the change of 1 • C in GMT t over this period, illustrating that climate change is more radically affecting extreme events than mean levels.
Next we consider summaries that also reflect the dependence structure of extreme events. There are no analytical closed form expressions of change. Instead we revert to using simulated fields of extreme events and presenting risk measures based on these. Our approach to simulation is to first simulate an event from the r-Pareto process, on Pareto margins, and then to back transform this to the data scale using the inverse of expression (9). As the r-Pareto process is stationary over time, we can generate identical events in Pareto margins to transform to any specific time using the time-varying Model 5 in expression (8). The r-Pareto process simulations are generated using the R package mvPot (de Fondeville et al., 2021), which exploits decomposition (13). Figure 8 shows simulated extreme events on the data scale under 2019 conditions and the exact same events in 1960 conditions (presented as a difference in value from 2019). By the nature of our modelling in Section 3.2, values of events below the threshold at a site do not change, but some events are changing by 3 • C almost everywhere.
Using such simulated events it is possible to generate Monte Carlo based estimates of risk measures. However it is inefficient to generate lots of small events when we are only interested in the largest ones. Here decomposition (13) proves hugely powerful as events can be scaled in size, through multiplying R with a known change in probability of occurrence. In particular, if we are only interested in events above temperature T • C for at least one site, then using the transformation (9)  We can determine a lower bound for R of b R = min s∈S v(s; T )/W max (s) for which we are sure not to miss generating any spatial field with a marginal value exceeding T • C. We therefore replace R by Rb R in decomposition (13), with R still Pareto generated, to gain efficiency. Despite this, some fields needed to be rejected if we are only interested in fields which somewhere in S exceed T • C. The simulation strategy set out above gives us replicated independent spatial fields which exceed T • C for at least one s ∈ S. We now propose risk metrics to summarise these. First consider a pairwise message of dependence, based on χ P but on the data scale, which combines the effect of changes in the marginal distributions over time with the estimated extremal dependence structure. Specifically, we define where s h is a randomly selected site in S with s h − s = h, i.e., the conditional probability of the observational process exceeding temperature T on day t at a site which is a distance h away from a site which does exceed temperature T on that day. Figure 9 presents estimates of this risk measure for T = 28−30 • C for each of 1960 and 2019 conditions. As expected extremal spatial dependence decreases with distance, but much more rapidly in 1960 than 2019. The effect of the critical temperature level T is particularly interesting as estimates of χ D o seem stable over T ∈ [28 − 30] for 2019, but decay with T over this temperature range for 1960. Thus, when viewed pairwise, the extent of spatial coverage of exceedances of critical levels is increasing over time, with that rate of change itself increasing with the value of the critical temperature level. Finally we look at a fully spatial risk measure in Figure 10. This figure plots the expected proportion of a temperature field on S, which exceeds T • C given that somewhere in the field T • C is exceeded, i.e., it plots E D o (T ; t) against T , defined as: where I(B) is the indicator function of event B. The figure shows the the expected spatial coverage of extreme events (as defined by a critical temperature occurring somewhere in Ireland) has increased from 1960 to 2019, with the relative change in spatial extent over this period increasing with temperature T . Even with events exceeding 28 • C, the spatial extent is more than doubling in this period.

Conclusion and Discussion
We have presented some new approaches to merge information from climate models into the spatial extreme value analysis of observed temperature. New methodological features include using outputs from an extreme value analysis of the climate model data to provide a covariate for the equivalent analysis of observational data, and dealing with r-Pareto processes in a missing data framework. We also presented novel metrics, combining both marginal and dependence features, to describe changes of spatial risk over time.
Our analysis was for Ireland daily maximum temperatures. We found that the climate model data were more helpful for marginal modelling of the observational data than for dependence modelling, as they have the potential to over-estimate extremal spatial dependence relative to the observational data. We pooled data from across stations to fit our model and find strong evidence that the covariates global mean temperature, Ireland's CO 2 emissions, and time were all needed to capture trends within current marginal data, but none of these appear to affect spatial extremal dependence (although we did not explore lagged relationships). We found that from 1960 to 2019 the occurrence rates of high threshold exceedances has doubled and extreme quantiles have increased by 1 − 3 • C, the latter 2 • C greater than global mean temperature. We also found that that events over thresholds that are critical for society have become much larger, approximately double, in spatial extent.
Section 5 illustrates the challenge of finding effective metrics to illustrate temporal change when considering extremes of spatial fields. For marginal variables such metrics, e.g., high quantiles, are well established and parsimonious. We see the development of spatial risk measures which enable the simple assessment of changing risk over time as an important avenue for further research. Spatial extreme value model inference also lacks well established diagnostics methods for assessing the fit with observed events. Pairwise, and potentially higher order versions of χ, used in Sections 4.2 and 5, can be helpful but are they sufficient? Winter et al. (2016) used severity-area frequency curves as a basis of comparison, but these focused only on assessing the performance of the model for the dependence structure.
Even with a suitable risk measure, when multiple covariates affect the risk then statistical methods are needed for determining the attribution of change to each covariate. Our analysis looks only at the change in extreme temperature events over the period of observations. Inference into the future is important, but that must allow for different scenarios their associated forecasting uncertainties. Our analysis used Irish CO 2 emission data, simply because it gave a better fit to observations than the equivalent global data, but that may not be a better predictor in the future.
Our methods for incorporating climate model data seem to be effective, but are they the most efficient way to extract the relevant information from climate model data? For example, could we learn about a spatially varying shape parameter for the observational data? Could the climate model data better inform us of temporal changes in spatial dependence? Could we learn more from different climate models or different ensemble members for a given climate model?
Although we assumed spatial stationarity in the analysis we present, use of the climate model data has helped identify non-stationarity, with the distance of sites from the coast also being an important descriptor of dependence. In exploratory work we have had reasonable success with the distance metric h ij between s i , s j ∈ S ofh ij = h ij + (δ 0 + δ 1 h ij ) exp{−δ 2 min(c i , c j )}, where h ij is the Euclidean distance between the sites, their respective shortest Euclidean distances to the coastline are c i and c j , and δ 0 , δ 1 and δ 2 are non-negative parameters. The distance metric can be incorporated in the Matérn covariance function (14). Future work would incorporate the resulting spatial non-stationary from using this metric. Our analysis focuses on the maximum daily temperature. Many other variables are important for assessing climate change for extreme meteorological events both marginally and jointly and over different time windows, e.g., rainfall over 1 hour and 24 hours. So the next step is to explore multivariate spatial extremes problems, again drawing knowledge from climate models to help improve extreme value inferences from observational data.