  Geostatistical models for multivariate applications such as heavy metal soil
contamination work under Gaussian assumptions and may result in underestimated
extreme values and misleading risk assessments (Marchant et al, 2011). A more
suitable framework to analyse extreme values is extreme value theory (EVT).
However, EVT relies on replications in time, which are generally not available
in geochemical datasets. Therefore, using EVT to map soil contamination
requires adaptation to be used in the usual single-replicate data framework of
soil surveys. We propose a bivariate spatial extreme mixture model to model the
body and tail of contaminant pairs, where the tails are described using a
stationary generalised Pareto distribution. We demonstrate the performance of
our model using a simulation study and through modelling bivariate soil
contamination in the Glasgow conurbation. Model results are given as maps of
predicted marginal concentrations and probabilities of joint exceedance of soil
guideline values. Marginal concentration maps show areas of elevated lead
levels along the Clyde River and elevated levels of chromium around the south
and southeast villages such as East Kilbride and Wishaw. The joint probability
maps show higher probabilities of joint exceedance to the south and southeast
of the city centre, following known legacy contamination regions in the Clyde
River basin.
