Exceedances and their uncertainty quantification

Greater awareness of human impact on the environment (e.g., IPCC, 2013) has led to an increase in the use of statistical modelling for predicting environmentally sensitive scientific processes of interest. For example, climate models help us to understand Earth's atmosphere, particularly those aspects associated with global climate change, such as the carbon cycle. Predictions are often made for the mean or median of a spatial process but, for risk assessment, strategic planning, and policy development, prediction of "extremes," or "exceedances" are also of interest; see Figure 1 for an example of rainfall-related exceedance probabilities over Australia.

Here, our research in environmental informatics is concerned with the problem of simultaneously predicting exceedances for a spatially dependent scientific process at many locations concurrently. This entails two challenges: Predicting the location, extent, and possibly the persistence of the exceedance, and then quantifying the prediction uncertainty; see here for a description of our approach to uncertainty quantification for hierarchical models. Further research in this area is given by Cressie and Suesse (2020).

Figure 1. The chance (percentage probability) of exceeding the median rainfall for October to December, 2015. The probabilities are generated from the Predictive Ocean Atmosphere Model for Australia (POAMA), the Bureau's dynamical climate model. Image from Australian Bureau of Meteorology, retrieved 23 Oct, 2015.

Definition of an exceedance

Consider a scientific process $Y(\bm{s};t)$, where $\bm{s}$ represents location and $t$ represents time. An exceedance region (a "hotspot" or a "coldspot") is defined as the set of locations $\bm{s}$ in a region $\mathcal{D}$ at time $t$, for which $Y(\bm{s};t)$ is above or below a given threshold, $u$.

Define a 'minimum' exceedance region ('coldspot'), $e^-(u)$, as the set of locations for which the spatio-temporal process $Y(\bm{s}; t)$ is lessthan $u$: \begin{equation} e^{-}(u) \equiv \{ \bm{s}\in D : Y(\bm{s}; t) < u \}; \end{equation} and define a 'maximum' exceedance region ('hotspot'), $e^+(u)$, as the set of locations for which the spatio-temporal process $Y(\bm{s}; t)$ is greater than $u$: \begin{equation} e^+(u) \equiv \{ \bm{s}\in D : Y(\bm{s}; t) > u \}. \end{equation}

Exceedances for CO$_2$ fluxes

One goal of climate models is to predict spatio-temporal changes in atmospheric carbon dioxide (CO$_2$) concentration, known as CO$_2$ flux $(\Delta_{CO2}(\bm{s}; t))$, and hence to predict regions of large positive CO$_2$ fluxes (sources), large negative CO$_2$ fluxes (sinks), and the movement of CO$_2$ between these regions. This is a challenging problem because sources and sinks are dynamic, and CO$_2$ flux is expensive and difficult to measure directly.

To predict CO$_2$ flux fields, climate models combine (1) predicted total-column atmospheric CO$_2$ concentrations (denoted $\hat{X}_{CO2}(\bm{s}; t)$), which are frequently obtained from remote sensing data (see Statistical remote sensing for a description of $\hat{X}_{CO2}(\bm{s}; t)$ prediction and our research in this area), with (2) the predicted creation, movement, and transformation of CO$_2$ in the atmosphere, which are obtained from highly complex, often deterministic chemical transport models such as Geos-Chem or LMDZ, and (3) a set of assumptions about the atmospheric state. See Regional flux inversion for a methane-based example of our research on predicting flux fields.

In our research, we use the resulting predicted flux field $\{\hat{\Delta}_{CO2}(\bm{s}; t): \bm{s} \in \mathcal{D}, t=1,2,...\}$ to predict the location, extent, and persistence of CO$_2$ exceedances in the atmosphere. For example, in Figure 2, a sink-type, or 'minimum,' exceedance region is defined as the locations for which the predicted CO$_2$ flux is less than $u=-6.59$ gC/m$^2$/day.

Figure 2. Predicted CO$_2$ flux field (gC/m$^2$/day) obtained from the inversion of a simulated prior flux field and atmospheric measurements (left panel); and the same predicted flux field with $e^-(u=-6.59)$ marked in green (right panel). Original flux predictions courtesy of F. Chevallier.

Exceedance predictions and their uncertainty

We wish to predict an exceedance for a spatio-temporal process $Y$, such as $\Delta_{CO2}(\bm{s}; t)$. However, due to measurement error and incomplete sampling, we cannot observe $Y$ directly. Instead we observe a prediction (or estimate) or a proxy for $Y$, such as a matrix of predicted flux fields obtained from satellite remote sensing data and a transport model. Define $\bm{Z} = (\hat{\Delta}_{CO2}(\bm{s}_i; t_j): i=1,\ldots,I; j=1,\ldots,J)'$ to be the predicted CO$_2$ flux field.

For a given level of statistical significance, $\alpha$, and a given exceedance threshold, $u$, we consider a $(1 - \alpha)$100% exceedance set for $e^{+}(u)$ or for $e^{-}(u)$ (henceforth, collectively referred to as $e(u)$), based on $\bm{Z}$, which we write as $e(\bm{Z},u; 1-\alpha)$.

To find the $(1 - \alpha)$100% 'minimum' exceedance set, we assume that all locations are NOT exceedances and, using a test statistic $T(\bm{s})$, we test the hypotheses, $H_0: Y(\bm{s}) \ge u$ versus $H_1: Y(\bm{s}) < u$, for each $\bm{s} \in \mathcal{D}$. The locations, $\bm{s}$, for which we reject $H_0$ are predicted exceedances, whilst other locations are not. In general, the predicted exceedance set $e(\bm{Z},u; 1-\alpha)$ is contained within the true exceedance set $e(u)$ with probability $(1-\alpha)$.

Methods for predicting exceedance regions must consider both the spatial dependence in the data and the multiple-testing problem. Prediction methods that account for spatial dependence include kriging-based approaches (e.g., Cressie, 1993), and hierarchical-model-based approaches (see here for a description of these). Alternatively, methods based on extreme value theory (Gumbel, 1958; Kharin and Zwiers, 2005; Cooley et al., 2007), max stable processes (Schlather and Tawn, 2003; Cooley et al., 2007), and Bayesian quantile regression (Jeffreys, 1961; Reich et al., 2011) are also used (see Davison et al., 2012, for a recent review of these methods).

The control of hypothesis-testing errors for large multivariate datasets is presently an active area of research. Historically, methods to control Type 1 error, such as the Bonferroni adjustment (Dunn,1961), were developed, and applied to spatial data (Adler, 1981). More recently, methods based on the False Discovery Rate (FDR) (Benjamini and Hochberg, 1995; Storey, 2003; Storey and Tibshirani, 2003), or the False Nondiscovery Rate (e.g., Craiu and Sun, 2008) have been developed and extended to correlated data (Benjamini and Yekutieli, 2001), and spatial data (Shen et al., 2002; Benjamini and Heller, 2007; Brown et al., 2014; and Sun et al., 2015). To account for multiple testing, recent approaches to prediction of exceedance regions have used a simulation-based approach (French and Sain, 2013), or a Bayesian approach (Bolin and Lindgren, 2015).

To obtain a $(1-\alpha)100\%$ predicted exceedance set, we use a simulation-based hypothesis-testing approach based on that used by French and Sain (2013). This approach applies fast, desktop computing to calculate the critical value used to reject $H_0$, and multiple testing at many locations is automatically accounted for.

Exceedance predictions for CO$_2$ fluxes

In our research, we have used predicted global flux fields from an observing system simulation experiment (OSSE), courtesy of Frederic Chevallier (Chevallier et al., 2005; Chevallier et al., 2007; Chevallier et al., 2010; Chevallier et al., 2014). Steps used in the OSSE to produce the predicted fluxes are given by Chevallier et al. (2007) as follows:

• Generate a set of pseudo-observations using a climatology of surface fluxes as boundary conditions to a transport model, such as LMDZ.
• Perturb the pseudo-observations consistently with the assumed error statistics.
• Perturb the surface flux climatology consistently with the assumed error statistics.
• Perform Bayesian inversion using the perturbed pseudo-observations as data and the perturbed climatology as the "prior field" to obtain a predicted flux field.
• Compare the predicted flux field to the known surface flux climatology to obtain error statistics.

From the OSSE, predicted flux fields for a global $97 \times 96$ grid (where the area of each grid cell depends on location), have been summarised for daytime (8am - 8pm local time) and nighttime (8pm - 8am local time), for four days per month during the period between January 2009, and January 2011. We fix a time $t$, and hence we drop the index $t$ from the notation. For the grid cell located at $\bm{s}$, a hierarchical-modelling approach was used to obtain a prediction, $\hat{Y}(\bm{s})$, for the hidden flux field, ${Y}(\bm{s})$, and the hypothesis tests described above were conducted using a test statistic of the form: \begin{equation*} T(\bm{s}) = \frac{\hat{Y}(\bm{s}) - u}{\sqrt{var(\hat{Y}(\bm{s}))}}. \end{equation*} Figure 3 gives an example of a 90% predicted exceedance set, for selected flux fields from the OSSE output, obtained using the simulation-based approach outlined above.

Figure 3. A 90% predicted exceedance set, $e(\bm{Z},u; 1-\alpha)$, based on the CO$_2$ flux field (gC/m$^2$/day) given in Figure 2, left panel. Notice that this is different from $e(u)$ shown in Figure 2, right panel, since it satisfies, $\text{Pr}(e(\bm{Z},u; 1-\alpha) \subset e(u))=(1-\alpha)100\%$. Original flux predictions courtesy of F. Chevallier.

Authored by S. Burden, 2015.