Joint Estimation of Extreme Precipitation at Different Spatial Scales through Mixture Modelling

by   Jordan Richards, et al.

Parsimonious and effective models for the extremes of precipitation aggregates that can capture their joint behaviour at different spatial resolutions must be built with knowledge of the underlying spatial process. Precipitation is driven by a mixture of processes acting at different scales and intensities. The specific process that drives the extremal behaviour of the aggregate will be dependent on the aggregate resolution; whilst high-intensity, spatially-localised convective events cause extreme high-resolution spatial aggregates, the contribution of low-intensity, large-scale fronts is likely to increase with the scale of the aggregate. Thus, to jointly model low- and high-resolution spatial aggregates, we require a model that can capture both convective and frontal events. We propose a novel spatial extreme values model which is a mixture of two components with different marginal and dependence models that are able to capture the extremal behaviour of convective and frontal rainfall. Modelling extremes of the frontal component raises new challenges due to it exhibiting strong long-range extremal spatial dependence. Our modelling approach is applied to fine-scale, high-dimensional, gridded precipitation data, where we show that accounting for the mixture structure improves the joint inference on extremes of spatial aggregates over multiple regions of different sizes.



page 10

page 14

page 29

page 30

page 31


Modelling Extremes of Spatial Aggregates of Precipitation using Conditional Methods

Inference on the extremal behaviour of spatial aggregates of precipitati...

Estimating high-resolution Red sea surface temperature hotspots, using a low-rank semiparametric spatial model

In this work, we estimate extreme sea surface temperature (SST) hotspots...

Identifying regions of inhomogeneities in spatial processes via an M-RA and mixture priors

Soils have been heralded as a hidden resource that can be leveraged to m...

Higher-dimensional spatial extremes via single-site conditioning

Currently available models for spatial extremes suffer either from infle...

A low-rank semiparametric Bayesian spatial model for estimating extreme Red Sea surface temperature hotspots

In this work, we focus on estimating sea surface temperature (SST) hotsp...

A Multi-Resolution Spatial Model for Large Datasets Based on the Skew-t Distribution

Large, non-Gaussian spatial datasets pose a considerable modeling challe...

Extreme event propagation using counterfactual theory and vine copulas

Understanding multivariate extreme events play a crucial role in managin...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

River flooding is caused by prolonged heavy rainfall over a designated catchment area, with the size of such areas varying dramatically between rivers in the UK. For example, the approximate catchment areas of the UK river Lune, that flows through Lancaster, and the Severn, the largest river in the UK, are and , respectively. Hence, a key aspect of river flood risk management is understanding the extremal behaviour of precipitation accumulations at different spatial scales. Whilst inference on precipitation extremes can be conducted at different spatial resolutions separately, and many studies have done so (see Institute of Hydrology (1999)

), ideally we require an approach that can be used for joint estimation of these variables at multiple spatial scales simultaneously. The benefits of such an approach include a framework for modelling dependence between aggregates over different regions and physically consistent inference, i.e., a natural ordering for quantiles of aggregates over nested regions

(Richards et al., 2021a). Moreover, this approach offers an increase in the reliability of inference on aggregates over lower-resolution spatial scales as more data will be used for modelling.

Building a parsimonious model for extreme precipitation at multiple spatial scales is difficult, as the important aspects of the underlying processes that drive the extremal behaviour of precipitation change with spatial resolution. Extreme precipitation is caused by a mixture of known and well-understood processes; convective storms are spatially localised events, over short time scales, and contribute to high intensity precipitation (Schroeer et al., 2018), whilst non-convective events, i.e., stratiform precipitation caused by frontal systems, are of much lower intensity but have a much larger area of effect and persistence (Berg et al., 2013, Catto and Pfahl, 2013, Gregersen et al., 2013). The distinction between the two types of events is an important consideration when modelling the upper-tail behaviour of spatial aggregates, particularly as the size of the region over which we aggregate increases; the contribution of frontal events is likely to increase with the size of the region.

Richards et al. (2021a) develop methodology for joint inference on precipitation extremes at different spatial scales. Specifically, they perform inference on the upper-tail of the aggregate variable


where denotes a spatial process for some spatial domain , and where is a sub-region of interest, with data available being realisations for times and sampling locations . To achieve this, they propose a model for the extremal behaviour of by adapting the spatial conditional extremes approach, first proposed by Wadsworth and Tawn (2019) and with extensions by Tawn et al. (2018), Shooter et al. (2019, 2021a, 2021b), Huser and Wadsworth (2020), Simpson and Wadsworth (2021), Simpson et al. (2021). Realisations taken from this model are used to conduct inference on for numerous regions simultaneously. Whilst their approach improved on previous approaches by ensuring self-consistent inference across different and flexible extremal dependence modelling, they make the restrictive assumption that both the marginal and dependence behaviour is stationary with respect to time; we know that for precipitation this assumption is poor. Here we adapt the modelling approach of Richards et al. (2021a); we propose separate models for convective and non-convective precipitation and simulate events from both models to approximate the upper-tails of for corresponding to different scales. Further extensions of the approach of Richards et al. (2021a) include an approach for simulating non-convective events, which suffer particularly badly from boundary effects due to their long-range dependence; a discussion of these complications and our solutions is given in Section 3.

We begin by assuming that times can be partitioned into two sets, denoted and , which correspond to “convective”, and “non-convective”, times respectively; that is, if , then the observed field is a convective event, and similarly for non-convective events. This classification is performed using an algorithm developed at the Met Office Hadley Centre, UK (see Section 4.2). We assume that there exist two different processes and , i.e., a convective, and non-convective, process respectively, and that these processes have different marginal and dependence structure. For each process, both the marginal behaviour and dependence structures are stationary with respect to the corresponding time sets; that is, the process is identically distributed to for all , and similarly for the non-convective process . We then define the mixture process by


for . As our interest lies in spatial, not temporal, aggregates, we thus drop the time index and let .

For and we adapt the approach proposed by Richards et al. (2021a), and in Section 2 propose slightly different parametric sub-models for the two processes, as we expect the mixture components to exhibit different extremal dependence structures. We quantify extremal dependence for a process using the upper tail index (Joe, 1997) for all , defined as


where denotes the distribution function of . The measure determines the extremal dependence class; if or , then and are asymptotically dependent or asymptotically independent, respectively (Coles, 2001). Richards et al. (2021a) find evidence that the process exhibits asymptotic independence at all spatial distances, i.e., for all , inferring that extreme events of occur increasingly more localised as the magnitude of the events gets larger. However, we find that exhibits asymptotic dependence at short-range and asymptotic independence at long-range.

Model inference is conducted by first estimating and for the observation times, which we do using the algorithm detailed in Section 4.2; we then fit separate models to the two classes of data using the inference framework described in Section 2.3.3. We simulate from our fitted models for and using a procedure described in Section 3, and then combine realisations from both models to estimate the upper-tail behaviour of . To illustrate the improvements to modelling that our approach provides to the approach detailed by Richards et al. (2021a), we apply both approaches to precipitation data from a convection-permitting climate model in Section 4. That is, we compare two modelling approaches for : an approach using the existing model, denoted throughout by , which ignores the mixture structure in ; and our improved approach where we use a mixture process with two components, i.e., , defined in (2). The two are equivalent when or .

2 Modelling extremes of each mixture component

2.1 Overview

Following Richards et al. (2021a), we adopt a two-step modelling approach. We first standardise the data using the marginal model described in Section 2.2; this model is presented for a generic stationary process , but is applicable separately for each of the processes and . The marginal fits are used to perform site-wise standardisation of data to standard Laplace margins; we denote these standardised processes by , and .

To characterise extremal dependence, we use the spatial conditional extremes framework proposed by Wadsworth and Tawn (2019); this is detailed for a generic standardised process in Section 2.3.1, but is applied, separately, to each of , and . Slight differences in modelling choices are made for the three processes, i.e., choices of parametric forms for some functions; discussion of these differences is provided in Section 2.3.2.

2.2 Marginal model

We model the site-wise marginals of using three components: GPD tails (Coles, 2001) above some high threshold , the quantile of , the empirical distribution for the bulk and a discrete mass at the lower tail, denoted , which corresponds to the probability that no rain occurs at site . The distribution function for , for all , is


where , and where , and and ; here denotes the empirical distribution function of strictly positive values of .

The marginal parameter functions and are represented through a basis of thin-plate splines. We estimate using a logistic generalised additive model (GAM) (Wood, 2006) and estimate for all using additive quantile regression (Fasiolo et al., 2021); the latter technique is computationally expensive and so we use only a subset of sites for estimation. Richards et al. (2021a) propose a technique for estimating by fitting a thin-plate spline through point-wise estimates of for each ; however, our approach provides more reliable estimates as more data are used for estimating the splines, rather than single quantile estimates. The scale parameter is estimated by fitting a GPD GAM to the exceedances of , i.e., , with a fixed shape parameter for all (Youngman, 2019). Our choice to fix over space was fully supported by our data and is a common approach taken when modelling spatial characteristics of extreme rainfall as it reduces the risk of having parameter identifability issues, see Thibaud et al. (2013), Zheng et al. (2015), Saunders et al. (2017) and Brown (2018).

Exploratory analysis of the data reveals that elevation is an important covariate in the marginal behaviour of precipitation; an observation supported by Coles and Tawn (1996), Cooley et al. (2007) and Cooley and Sain (2010). Hence, we allow and to vary with both location and elevation; separate spline bases are used for both. Each spline consists of as few knots as possible, i.e., four, to ensure that the splines are smooth; this is to make the marginal fits physically interpretable and to reduce over-fitting. Our approach for inference differs from that proposed by Richards et al. (2021a) as we allow the margins to change with elevation.

2.3 Dependence model

2.3.1 Generic model

We model extremal dependence in the standardised process , given that it is extreme for some , by conditioning on the process being above some high threshold at a specified site for each . The process is assumed to be stationary with dependence a function of distance, i.e., for , and some distance metric . Following Wadsworth and Tawn (2019), we assume that there exists normalising functions , with , and , such that for each , as ,


where is a standard exponential variable, the process is non-degenerate for all where and it is independent of . That is, there is convergence in distribution of the normalised process to , termed the residual process and almost surely.

In terms of a parametric specification of , we also follow Wadsworth and Tawn (2019) and let


where and . This function determines the strength and class of extremal dependence within , see (3); allowing asymptotic dependence up to distance from , and asymptotic independence thereafter, with strength of dependence decreasing with . Setting and for all makes an Pareto process, which have been used for aggregated precipitation modelling by de Fondeville and Davison (2020).

Richards et al. (2021a) impose that the residual process for has delta-Laplace margins with location, scale and shape parameters determined by the spatial functions , and , respectively. We denote this by , where the delta-Laplace density is for and and , where denotes the standard gamma function and . Dependence within , for each , is modelled using a standard form Gaussian process with a stationary Matérn correlation function, denoted . We adopt the same parametric forms for , , , and as Richards et al. (2021a), unless stated otherwise in Section 2.3.2. These parametric forms allow the extremal dependence model for to exhibit independence at long-range as , , , and when .

The function accounts for potential anisotropy in the extremal dependence structure of . As Richards et al. (2021a), we define the distance metric , which involves an elliptical transformation of and then finds the great-circle, or spherical, distance between them; the metric is parametrised by and which control the rotation and coordinate stretching effect, respectively, with corresponding to isotropy.

2.3.2 Differences for non-convective component

We apply the model described in Section 2.3.1 to each of , and . For and , we use exactly the same parametric forms for the dependence functions as proposed by Richards et al. (2021a); for , we adapt the and functions, denoting these new forms and .

Richards et al. (2021a) advocate the use of a function for , say , that satisfies the property . They illustrate that this, combined with small values of for small , allows the model to exhibit short-range spatial roughness; a desirable property for a process generating convective precipitation events. As we expect to generate much smoother spatial events, we do not require the property that ; to support this claim, we illustrate observations of and in Figure 2. We instead adopt the approach of Shooter et al. (2021b) and let


so for all and as with . In our application, we consider both forms of for , but found that (7) provided a better fitting model.

The scaling function gives long-range independence. However, non-convective events can have a very large spatial extent, e.g., up to (Houze Jr, 1997). Long-range independence is not required for when the domain of interest is small enough, as in our application in Section 4. To ensure a better fitting model for , we let the function be


which is equivalent to when .

2.3.3 Inference

Inference for each of the three extremal dependence models is conducted using the pseudo-likelihood procedure proposed by Richards et al. (2021a), which they illustrate works well for models of this type and for application to precipitation data. This technique requires selection of a sub-sample of triples of sampling locations, subject to the maximum pairwise distances subceeding a value of . Maximisation of a pseudo-likelihood is then conducted using observations of the process at these triples via a triple-wise censored likelihood. We require a censored likelihood to handle point masses in the marginals of caused by zeroes in the marginals of . The censoring threshold for , we denote , is found by transforming to the Laplace scale, i.e., , where denotes the standard Laplace distribution. We apply this approach and use different values of and for each of the three processes; further details are provided in Section 4.4. Inference is conducted for model (5) by fixing a high threshold and assuming that the limit holds in equality; the value of differs between the three processes. Finding a threshold such that and , defined in (5), are independent may be infeasible with mixtures of different dependence present in , e.g., for . This provides further support for our use of separate extremal dependence models for the mixture components, as the assumptions made in (5) are more likely to hold when applied to and , rather than , and at lower thresholds ; we find the latter to be the case in our analysis.

3 Simulating events

3.1 Efficient simulation

To simulate from the fitted models for , and , we need to extend the procedure proposed by Richards et al. (2021a) due to long-range dependence features of frontal rain; we present this technique for a generic and , but this approach can be applied to any of the three considered processes. The approach we propose is tailored to a specific aggregate region , but our interest may lie in aggregates over multiple regions for ; to ensure that we can perform joint inference on aggregates at multiple scales, we let throughout.

Richards et al. (2021a) discuss edge effects that occur when considering aggregates over regions , as the simulation procedure never generates events for which the conditioning sites lies outside of the boundaries of but an extreme event is observed within . To address this issue, they proposed constraints for choosing ; they require that the boundaries of are at least distance within the boundaries of . They choose the value of so that there was a predetermined small probability of observing a large event at any given that there is an extreme event at any single site outside of . When exhibits particularly strong extremal spatial dependence, picking a suitable may not be feasible; in our application we find this to be the case for . To address this issue, we introduce a set of conditioning sites that satisfy and create realisations of for . Whilst Richards et al. (2021a) set , here, when it is possible to sample from locations outside of , we allow for cases where , then using we observe a reduction in edge effects when are sufficiently far from .

We now assume that is a set of dense discrete sites, i.e., the set of observed sampling locations . When is large, it may be computationally infeasible to simulate fields for all , then instead we can simulate fields on a sufficiently dense sub-region , where . Illustrations of , , and for our application are presented in Figure 1

and a heuristic for choosing both

and is given in Section 3.2; these sets vary with the process considered.

Figure 1: Left: A map of elevation for the spatial domain of interest. Right: Regions , and . Aggregate regions with corresponding areas are coloured red, green, blue, cyan; regions include both the coloured and interior points and are numbered 1 to 4 in Figures S8 and 6. The orange and black points denote and , respectively; the purple points outside of the boundaries of denote . Note that , , and ; these values mean that and .

To simulate , we draw realisations of


for , with probability


and otherwise draw realisations of


from the observed . To draw realisations of (3.1), we use Algorithm 1.

  1. For with :

    1. Draw a conditioning location from with uniform probability .

    2. Simulate and set .

    3. Simulate a field from the residual process model.

    4. Set for each .

  2. Assign each simulated field an importance weight of


    for , and sub-sample realisations from the collection with probabilities proportional to these weights.

  3. Transform each to using the marginal transformation (4). If , set , where for some , is above its -th quantile.

Algorithm 1 Simulating conditional spatial fields corresponding to formulation (3.1)

Using Algorithm 1, we can draw realisations of and ; we can then use these to derive samples of by drawing from the samples of and with probabilities and , respectively, see (2). Estimates of and are both derived empirically; the former from the data and the latter from the simulated replicates. We denote estimates of derived using and by and , respectively.

3.2 Choosing and

Although we have taken to be a discrete set, the heuristic we describe for choosing and can be extended to a spatially continuous setting. We begin by considering with and then build dependent on . For Algorithm 1 Step 12, we require the approximation


to account for the conditioning event in (3.1); we can improve the accuracy of this approximation by ensuring that sites in are sufficiently spread out across . We first set as we expect extreme events within this set to have the largest effect on the tail behaviour of . If this choice of does not give a good approximation (13), then we additionally sample sites uniformly at random across and add these to .

In our application, we found that setting was a reasonable choice to make for and , as no improvement in our inference for was achieved by using . However, this was not true for , which exhibits strong extremal dependence at even the greatest distances within . For such processes, we need to increase to ensure that we have conditioning sites that are sufficiently far from such that simulated do not contribute to the extremes of . We take and let , with a set of “fake” sampling locations without any observations and . We create via a brute-force approach; we add independent Gaussian noise to the coordinates of a site in , rejecting new points .

In both cases, and should be chosen as large as is computationally feasible. Suitable values for the constants and

, and the variance of the noise term used to create

, can be chosen through validation techniques, such as the simulated aggregate diagnostics (see Figures S8 and 6). If such fits look poor, one can increase the value of these constants to improve the fits.

4 Application

4.1 Data

Similar to that of Richards et al. (2021a), we consider average hourly precipitation rates (/hour) taken from the convection-permitting climate model UKCP18 (Lowe et al., 2018), 1980-2000. However, the data here are on the finer resolution native climate model grid over a much larger domain and the dimension is increased from to . We analysed the first and fourth ensemble members from the UKCP18; however, we present the latter analysis only, as we identified abnormally large and physically unrealistic values in the former. Each observation corresponds to the average over the assigned spatio-temporal grid-box. The sampling locations are grid boxes and the spatial domain of interest is the region of the UK, approximately centred at Northampton and covering an area of , pictured in Figure 1, with data only sampled over land. As extreme hourly precipitation is more likely to occur in summer, we use June-August observations only, fields, to avoid having to handle seasonality effects. We follow Richards et al. (2021a) and treat the centre of each grid box as a sampling location, and use the great-circle distance as our distance metric; we also set all values of the data less than /hour to zero111This level of precipitation would be recorded as zero by a rain gauge..

4.2 Mixture component classification

To classify observation times as convective or not for the observed fields

, we use Algorithm 2 developed at the Met Office Hadley Centre, UK (Kendon et al., 2012, Roberts and Kendon, 2020) for summer precipitation data with gridded sampling locations. The algorithm classifies the observed field at time as convective by using the gradient of the surface of in a neighborhood around for all . This approach is a good classifier for convective rainfall as these events are highly localised, so large gradients are expected if convection is present. Gradients are calculated as the difference between precipitation rates at adjacent grid-boxes in a grid, denoted , centered at . If it is not possible to create , e.g., at the edges of , then is removed. Any field not identified as convective is automatically classified as non-convective; this includes any fields with no rain at any . Algorithm 2 involves four hyper-parameters: and , we take to be and and ; these were tuned by the Met Office specifically for the climate model we use.

For all :

  1. For all :

    1. Identify the local neighbourhood for , defined in Section 4.2.

    2. Evaluate all differences .

    3. Calculate the proportion .

    4. If , then is labelled as convective and hence .

  2. If none of for all , are labelled as convective, then .

Algorithm 2 Identify convective fields

Algorithm 2 gives and , so , with examples of classified fields of and presented in Figure 2. These selected fields were randomly sampled from fields which gave componentwise maxima for their respective processes. Fields identified as being non-convective appear smoother over space and with much lower marginal magnitude than the convective fields; note the difference in the scales of Figure 2.

Figure 2: Observed extreme fields identified as convective (top) and non-convective (bottom) (/hr).

4.3 Marginal analysis

Marginal analysis is conducted by fitting the model described in Section 2.2 to each of the three datasets, i.e., convective, non-convective, and all. We set in (4) for all three processes and use a subset of sites sampled randomly over to estimate ; whereas to estimate , and , we use all sampling locations. Figures S1-S3 in the Supplementary Material B (Richards et al., 2021b) give estimates for the parameters of the marginal models and -year marginal return level estimates for , and , respectively, with the latter being estimated for the specified process only. Different image colour scales are used across each panel and each of Figures S1-S3. We observe similar patterns in estimates of and for each of the three processes, with and decreasing and increasing with elevation, respectively.

There are differences between the estimates of and the 20-year return levels for the three processes; for and , we observe spatially smooth estimates of both, with larger values being found in the east of the domain (see Figures S1 and Figures S3) for both increase with elevation, indicating that the physical processes driving convection are only weakly affected by orthography, unlike frontal rain (see Figure S2). We observe much higher -year return levels for and than for , which is consistent with their different physical properties. This property also holds for all larger return levels since the shape parameter estimates for , and , are 0.226, -0.075 and 0.287, respectively.

To assess the fit of the GPD GAM models, for each of , and we present Q-Q plots of the marginal fits at five randomly sampled locations in Figures S4-S6, Appendix B, respectively. All figures show good individual fits for each of the processes. To evaluate the fit over all locations, we use a pooled Q-Q plot (Heffernan and Tawn, 2001)

, transforming all data onto standard exponential margins using the fitted model; these plots are also given in these figures. Again, we observe good fits for each process. Confidence intervals for the Q-Q plots are estimated using the following bootstrap procedure: we create

bootstrap samples of the data using the stationary bootstrap approach of Politis and Romano (1994) with expected block size of 48 hours. With treated as fixed across all samples, we then estimate and for each bootstrap sample. For the pooled diagnostic plot, we apply the marginal transformation to the original data using the estimated GPD parameter sets.

4.4 Dependence modelling

We proceed by fitting the extremal dependence models described in Section 2.3 separately to , and . We use a different exceedance threshold, i.e., in (5), for each process; the quantile for and the quantile for and . The thresholds were chosen so that the limiting assumptions in (5) were applicable; other choices were considered but these led to poorer inference on the extremal behaviour of spatial aggregates, which we quantified using the measures and diagnostics described in Section 4.5. Ideally, for we would set such that the number of observations used for inference was identical to the mixture approach to provide a fair comparison. However, we found that using such a threshold led to poorer inference for , with significant evidence that assumptions in (5) fail for any lower threshold choice.

Inference is conducted using the stratified sampling regime described in Section 2.3. We use the same value of for all three fitted models; we follow Richards et al. (2021a) and set . Different values of were also taken: for and and for , with the latter reflecting the long-range extremal dependence of . Inference for all three extremal dependence models was conducted using different sub-samples of sampling locations. For we tried both forms of , i.e., and in (7); with giving a better fit, we present findings for this form only. Parameter estimates for all models are given in Supplementary Material A (Richards et al., 2021b). For and , we found that fixing and did not restrict the quality of either fit.

In Figure 3, we fix a conditioning site in the centre of the domain and evaluate all dependence functions for each process at each pairwise anisotropic distance, i.e., for and ; these are then plotted against the great-circle distances in the original coordinate system, denoted . Figure 3 indicates that there are similarities within the extremal dependence structures of and , as most of the functions are approximately equal with slight differences for and . In contrast, there are widely different estimated structures for the and functions for , with decaying much slower and independence not being achieved for any pairs of sites within .

The estimate of given in Table S1 suggests that the process is asymptotically dependent up to a distance of in the anisotropic setting, corresponding to roughly in the original coordinate system. Evidence of asymptotic dependence, even at short-range, has not been found previously using the spatial conditional extremes modelling approach, see e.g., Wadsworth and Tawn (2019), Huser and Wadsworth (2020), Simpson and Wadsworth (2021) and Shooter et al. (2021a). As extremal dependence is predominantly exhibited through the function, this suggests that extreme realisations of are much smoother than extreme events drawn from the other two processes, which agrees with our understanding of convective and non-convective precipitation. Keef et al. (2013) identify additional constraints for the Heffernan and Tawn (2004) framework, which in the spatial context correspond to whenever ; our fitted model breaks this constraint with for , but this minor deviation is not problematic.

Figure 3: Estimates of extremal dependence functions evaluated at for , i.e., anisotropic distances, against great-circle distances , which are given in . The conditioning site is in the centre of the spatial domain . The colours correspond to the estimates for the different spatial processes; these are green, red and blue for and , respectively.

To compare the model fits, i.e., and , we propose the following diagnostic that we present for a generic process. We investigate how the “conditional process” changes with distance for , where denotes the -year return level for . Figure 4 gives point-wise estimates for the median, and the and marginal quantiles, for and for the conditional process; these are obtained from realisations from the fitted model. To achieve the largest possible distances, is chosen to be on the boundary of and to ease comparison, we consider a transect of points in only. We take into account the respective length of observation periods when evaluating the respective year return levels, e.g., for we take and similarly for and . The diagnostics for and are almost identical, with and approximately equal, and with the conditional medians of the respective processes both decaying quickly with distance. For , the conditional median decays fairly slowly, with non-zero values at even the largest distances.

Figure 4: Summary statistics for against distance with and in the left and right plots, respectively; the process is considered along a transect of points only. Solid lines correspond to estimates for conditional medians, dashed lines denote confidence intervals. Lines are coloured green, red and blue for , and , respectively. Centre: red and blue points denote the transect and , respectively.

Figure 5 presents realisations of for and , where and each was uniformly sampled over , with their locations identified; similar plots are given for in Figure S7, Supplementary Materials (Richards et al., 2021b). There are similarities for and , with characteristics we anticipated for convective rainfall, i.e., high intensity, spatially localised events with a large proportion of the domain being dry. In contrast, produces events that are much smoother spatially, cover a much large area, and are lower in their intensity.

Figure 5: Extreme precipitation fields (mm/hr). Realisations from the fitted models for , with the top and bottom rows corresponding to and , respectively. The conditioning sites are given by the red crosses. Scales differ within each panel and row.

4.5 Inference on spatial aggregates

For each process, we draw realisations using Algorithm 1, being increased until the estimated quantiles of were stable. For and , we used ; for , . Regions and , described in Section 3.2, are illustrated in Figure 1; and are not changed over the processes, but we take for and and for we create using the heuristic described in Section 3.2, with illustrated by the purple points in Figure 1. Using these regions, we create samples of for different regions (see Figure 1) by drawing from and with estimated probability and , respectively.

Recall that denote samples from created using our new modelling approach and denote samples using the single process approach of Richards et al. (2021a). To compare mixture and non-mixture approaches, we present Q-Q plots in Figures 6 and S8 (Supplementary Materials (Richards et al., 2021b)) for spatial aggregates. In Figure S8, we assess the fit of individual mixture components and by using Q-Q plots of the variables