Coexchangeable Process Modeling for Uncertainty Quantification in Joint Climate Reconstruction

Abstract

Any experiment with climate models relies on a potentially large set of spatio-temporal boundary conditions. These can represent both the initial state of the system and/or forcings driving the model output throughout the experiment. These boundary conditions are typically fixed using available reconstructions in climate modeling studies; however, in reality they are highly uncertain, that uncertainty is unquantified, and the effect on the output of the experiment can be considerable. We develop efficient quantification of these uncertainties that combines relevant data from multiple models and observations. Starting from the coexchangeability model, we develop a coexchangeable process model to capture multiple correlated spatio-temporal fields of variables. We demonstrate that further exchangeability judgments over the parameters within this representation lead to a Bayes linear analogy of a hierarchical model. We use the framework to provide a joint reconstruction of sea-surface temperature and sea-ice concentration boundary conditions at the last glacial maximum (23–19 kya) and use it to force an ensemble of ice-sheet simulations using the FAMOUS-Ice coupled atmosphere and ice-sheet model. We demonstrate that existing boundary conditions typically used in these experiments are implausible given our uncertainties and demonstrate the impact of using more plausible boundary conditions on ice-sheet simulation. Supplementary materials for this article are available online, including a standardized description of the materials available for reproducing the work.

Keywords:

• Bayes linear methods
• Exchangeability analysis
• Multi model ensemble

1 Introduction

Numerical experiments are vital tools to climate science. Knowledge of physical processes embedded into software can simulate events that are otherwise impossible to observe at the required spatial and temporal resolutions. Most climate simulators use boundary conditions to represent the non-computed physical processes on which the simulator relies. Generally, boundary conditions are fixed using a reference run or runs from existing multi-model ensembles (MMEs) that explicitly model the boundary condition’s physical process. For example, to run an ice-sheet model over the last glacial maximum (LGM) requires climate variables such as temperature, and precipitation (Gregoire, Payne, and Valdes 2012) that can be obtained using the Paleoclimate Model Intercomparison Project (PMIP) ensemble of simulator runs (Ivanovic et al. 2016). MMEs can be biased and do not necessarily span the uncertainty of the physical process (Salter et al. 2022).

If modeling boundary conditions is viewed as a statistical reconstruction problem, there is a rich literature in statistics that attempts to combine data from multiple models and historical observations to infer spatio-temporal climate properties. Rougier, Goldstein, and House (2013) present a generalized framework, therein termed the coexchangeable model, where exchangeability judgments over an MME along with an assumption that each model in the ensemble has the same a priori covariance with the field they aim to represent lead to a simple model that can be used to estimate a true unobserved process. Second order inference via Bayes linear methods (Goldstein and Wooff 2007) for this model ensures scalable reconstructions for a spatio-temporal field, and the requirement for only prior means and variances implies an easier prior modeling task. The coexchangeable model has seen some interest in more recent research, such as in Sansom, Stephenson, and Bracegirdle (2021) where it is used to model emergent constraints for future climate projections. Alternative to the methodology of Rougier, Goldstein, and House (2013), a line of research stemming from Chandler (2013) proposes a similar fully probabilistic model with individual weightings of each simulator rather than via an assumption of exchangeability. As noted in Rougier, Goldstein, and House (2013), the model in Chandler (2013) is a special case of the coexchangeable model. Other statistical reconstructions of boundary conditions exist; however, they are either derived from simulator data alone and are limited by the span of the ensemble (e.g., Salter et al. 2022), or are reconstructed from large and dense datasets (e.g., Liu and Guillas 2017; Zhang and Cressie 2020) and so are inappropriate here given the relative data sparsity at the LGM. There are a family of post-processing methodologies in the climate sciences, therein termed offline data assimilation, that adjust model output with direct observations or proxy data based on ensemble methods (e.g., Steiger et al. 2014). We show that these methodologies, as well, may be subsumed by the coexchangeable model. Online data assimilation techniques are not possible for large climate ensembles: the simulators are run on large supercomputers, are highly bespoke, and very limited access to the model architectures is given outside of the modeling groups.

Joint reconstruction of two or more physically distinct fields may, from a theoretical perspective, appear straightforward within hierarchical Bayesian frameworks such as those proposed by Chandler (2013). However, scalability becomes a much bigger challenge. This is particularly so when the fields are highly dependent, as many in climate are. In this study we jointly model sea-surface temperature (SST) and sea-ice concentration (SIC) to drive a coupled atmosphere and ice-sheet model, where the dependence between the boundary conditions is critically important. The need for scalable inference and generality makes a coexchangeable approach desirable; though extensions are required to incorporate the additional assumptions of conditional exchangeability between the fields. Thus, we develop a coexchangeable process model that offers a Bayes linear analogy to the natural Bayesian hierarchical model for the problem. From simple and natural exchangeability judgments, we develop efficient, scalable inference for joint reconstruction.

This article proceeds as follows. Section 2 provides background and context to the applied modeling problem. Section 3 reviews the existing statistical literature on Bayes linear statistics and exchangeability analysis. Section 4 extends the coexchangeable model for coupled processes, providing a Bayes linear analogy to the Bayesian hierarchical model for which we then present scalable inference via geometric updating. Section 5 uses the developed methodology to reconstruct SST and SIC. Section 6 discusses the results of using the reconstructions of SST and SIC to simulate ice sheet-atmosphere interactions at the LGM, and Section 7 concludes. Code and data are provided as part of an R package downloadable from github.com/astfalckl/exanalysis.

2 Reconstructing Boundary Conditions of SST and SIC

Modeling and understanding paleoclimate events is crucial in understanding potential effects of future climate change: one way to calibrate climate and ice-sheet behavior is to simulate the past (Kageyama et al. 2017; Harrison et al. 2015). The last major deglaciation (∼21–7 kya) is the most natural period to study, given the relatively rich source of observations on the climate and ice sheets compared to the more distant past. To begin to study deglaciation, an ice sheet needs to be grown within the model under the steady state boundary conditions. The ice-sheet is grown by forcing an ice-sheet model with the relatively stable conditions of the LGM (∼23–19 kya) until convergence. Seasonal variation is known to play a significant role in this process (Joughin, Alley, and Holland 2012), so the boundary conditions contain the seasonal cycles. The character of the resulting ice sheets are sensitive to these boundary conditions and so it is crucial to use an accurate reconstruction, realistic uncertainty quantification and a method for perturbing the boundary conditions under uncertainty to force the ice sheet model.

Existing reconstructions of LGM SST and SIC are predominantly derived from either paleodata syntheses or via numerical simulation. Data only reconstructions include CLIMAP (CLIMAP Project 1981), GLAMAP (Sarnthein et al. 2003), MARGO (Kucera et al. 2005), and more recently Paul et al. (2020). As these global reconstructions are based solely on proxy-based paleodata, they are subject to large measurement error, biases in polar regions, incomplete spatial coverage, and poor temporal resolution. Coupled simulations of the LGM are useful to ensure spatio-temporal coverage and consistency of SST and SIC, but these can also be very different from observations in critical regions for growing ice sheets (Salter et al. 2022). PMIP is the most notable experimental body that guide protocols for coupled ocean-atmosphere models to simulate the LGM climate (Kageyama et al. 2017, 2021). The PMIP community has produced MMEs for different phases of the project run with different generations of model and slight adjustments in inputs. It is common practice to use a PMIP simulator output to directly force ice sheet simulations (e.g., Gregoire et al. 2016). One study to use both paleodata and simulators is Tierney et al. (2020) who adjust a model with observations to provide SST reconstruction, with uncertainty, based on a single simulator. Differences in model physics typically induce more variability than perturbations to the parameters of a single model. Therefore, reconstructions based on a single model may be biased and overconfident. Our methodology allows for the first joint reconstruction of SST and SIC that coherently combines the PMIP models (we use iterations PMIP3 and PMIP4) with available proxy data to deliver boundary conditions with uncertainty quantification.

We use three sources of data to inform our reconstructions: PMIP simulations, SST proxy data, and maximum sea-ice extents. As with Rougier, Goldstein, and House (2013) and Sansom, Stephenson, and Bracegirdle (2021), we select a single representative simulation from each modeling group that contributed to the PMIP3 and PMIP4 MMEs in order to make the assumption of prior exchangeability reasonable. We use the MARGO SST proxy data compilation where sea-bed sediment core samples were used to infer the true LGM SST, supplemented with some more recent data from Benz et al. (2016). Note, many of the SST measurements are far from any icesheets but can still useful in constraining a global climate simulation. Finally, we use the Southern Hemisphere maximum sea-ice extent as published in Gersonde et al. (2005); Northern Hemisphere extents are only available for specific regions (e.g., De Vernal et al. 2005; Xiao et al. 2015), and so we use a simple estimate of the Northern Hemisphere sea-ice extent provided by subject matter experts. SST proxy data and maximum sea-ice extents are shown in Figure 2 in Section 5 alongside the specification of the statistical model.

3 Existing Statistical Methodology

SST and SIC are dependent quantities: warmer waters will support less sea-ice and vice-versa. The physics that govern the relationship between SST and SIC is represented by a series of partial differential equations, the structure and parametrizations of which change across the simulators and with reality. For example, certain simulators are predisposed to supporting more or less sea-ice at a given SST. The dependencies between SST and SIC are naturally modeled as an emergent property of the model and captured within a hierarchical framework using conditional probability statements. However, in climate modeling, specification of the probability distributions is not obvious and computation for large climate models is prohibitive. An alternative view treats expectation, rather than probability, as the primitive quantity (De Finetti 1975); probabilities are then the expectations of indicator functions for events. This motivates second order approaches such as Bayes linear methods (Goldstein and Wooff 2007), where inference concerns expectations and variances directly rather than as a by-product of probabilistic inference. Model specification thus only concerns specifying expectation and variance, rather than the full probability distributions. Rougier, Goldstein, and House (2013) show the advantages of second order specification for climate modeling and further show that based on certain judgments of exchangeability, efficient methods for inference of MMEs may be formed. Our application requires extending the theory of Rougier, Goldstein, and House (2013) to exchangeable processes, that is, second order hierarchical models. First, we review the existing methodology: Section 3.1 provides a short introduction to Bayes linear theory, Section 3.2 defines second order exchangeability, and Section 3.3 presents the coexchangeable model of Rougier, Goldstein, and House (2013).

3.1 Bayes Linear Theory

Under the Bayes linear paradigm, beliefs on random quantities are described via expectations and variances and are then adjusted by data. The belief specifications define an inner product space in which the random quantities live; the inner product space is analogous to a probability distribution in probabilistic Bayesian analysis. Consider random quantity, $\mathrm{X}$, with observations ${\mathrm{X}}^i$, defined on the Hilbert space $\mathcal{X}$ endowed with inner product $〈\mathrm{X},\mathrm{Y}〉=\mathbb{E}[{\mathrm{X}}^{⊺}\mathrm{Y}]$. Denote by $\mathrm{D}$ the collection of observed data, $\mathrm{D}=({\mathrm{X}}^1,\dots ,{\mathrm{X}}^m)$, as a concatenated vector of the observations. In general, the ${\mathrm{X}}^i$ do not necessarily have the same length or a-priori belief specifications and can represent multiple sources of data that inform $\mathrm{X}$. The random quantity of interest, $\mathrm{X}$, may be multivariate or univariate, in which case its inner product is simply $\mathbb{E}[\text{XY}]$. Later, $\mathrm{X}$ will denote the unknown spatio-temporal field of SST, and the ${\mathrm{X}}^i$, the observed simulator outputs over which we will assume exchangeability. We write the adjusted expectation in terms of $\mathrm{X}$ and $\mathrm{D}$ as ${\mathbb{E}}_{\mathrm{D}}[\mathrm{X}]$, that is, the expectation of our beliefs $\mathrm{X}$ adjusted by data $\mathrm{D}$. Adjusted expectation is defined as the element in the subspace $\mathcal{D}=\text{span}[1,\mathrm{D}]$ that minimizes $\left|\right|\mathrm{X}-{\mathbb{E}}_{\mathrm{D}}[\mathrm{X}]\left|\right|$ and has solution (1) $${\mathbb{E}}_{\mathrm{D}}\left[\mathrm{X}\right]=\mathbb{E}\left[\mathrm{X}\right]+\text{cov}[\mathrm{X},\mathrm{D}]\text{var}{\left[\mathrm{D}\right]}^{†}(\mathrm{D}-\mathbb{E}[\mathrm{D}])$$(1) where $\text{var}{[\mathrm{D}]}^{†}$ is any pseudo-inverse of $\text{var}[\mathrm{D}]$, most commonly the Moore-Penrose inverse. Equation (1)(1) $${\mathbb{E}}_{\mathrm{D}}\left[\mathrm{X}\right]=\mathbb{E}\left[\mathrm{X}\right]+\text{cov}[\mathrm{X},\mathrm{D}]\text{var}{\left[\mathrm{D}\right]}^{†}(\mathrm{D}-\mathbb{E}[\mathrm{D}])$$(1) describes the orthogonal projection of each element of our beliefs, $\mathrm{X}$, onto $\mathcal{D}$. The adjusted variance, ${\text{var}}_{\mathrm{D}}[\mathrm{X}]$, is the outer product $\mathbb{E}\left[(\mathrm{X}-{\mathbb{E}}_{\mathrm{D}}[\mathrm{X}]){(\mathrm{X}-{\mathbb{E}}_{\mathrm{D}}[\mathrm{X}])}^{⊺}\right]$ given ${\mathbb{E}}_{\mathrm{D}}[\mathrm{X}]$ in (1) and is (2) $${\text{var}}_{\mathrm{D}}\left[\mathrm{X}\right]=\text{var}\left[\mathrm{X}\right]-\text{cov}[\mathrm{X},\mathrm{D}]\text{var}{\left[\mathrm{D}\right]}^{†}\text{cov}[\mathrm{D},\mathrm{X}].$$(2)

Derivations of (1) and (2) are found in secs. 12.4–12.5 of Goldstein and Wooff (2007).

3.2 Bayes Linear Analysis of Exchangeable Data

Exchangeability for a sequence of random quantities within a probabilistic Bayesian analysis represents a simple a priori symmetry judgment that implies that any finite sub-collection of quantities within the sequence have the same distribution. Second order exchangeability for such a sequence is an analogue for judgments when expectation is primitive, and implies that any finite sub-collection of quantities share the same joint inner product space, that is prior expectation and variance. If data, $\mathrm{D}$, are second order exchangeable we may, according to the second order representation theorem (Goldstein 1986), write (3) $${\mathrm{X}}^i=\mathcal{M}+{\mathcal{R}}^i$$(3) for a common mean $\mathcal{M}$ and uncorrelated residuals ${\mathcal{R}}^i$. Second order exchangeability implies that all observed ${\mathrm{X}}^i$ are of the same length, $\text{cov}[{\mathrm{X}}^i,{\mathrm{X}}^{i\prime }]=\text{cov}[\mathrm{X},\mathrm{X}\prime ]$ and $\text{var}[{\mathrm{X}}^i]=\text{var}[\mathrm{X}],\forall i,i\prime$. For second order exchangeable sequences there is predictive sufficiency for updating beliefs on $\mathrm{X}$ by only updating $\mathcal{M}$ by the data $\mathrm{D}$. Geometrically, this means that given $\mathcal{M},\mathrm{D}$ and $\mathrm{X}$ are orthogonal and thus uncorrelated; these results are established in Goldstein (1986). Further, the sample mean, $\overline{\mathrm{X}}=\frac{1}{m}{\sum }_{i=1}^m{\mathrm{X}}^i$, is Bayes linear sufficient for updating beliefs on $\mathcal{M}$, and consequently on $\mathrm{X}$. Second order exchangeability affords two main advantages: first, belief specifications are simplified; and second, sufficiency of the sample mean makes inference independent of the number of samples m.

3.3 Exchangeability Analysis for Multi-Model Ensembles

Rougier, Goldstein, and House (2013) leverage second order exchangeability to describe a Bayes linear approach for modeling MME’s. In what follows, we explicitly define multivariate quantities, represented by the bold font. Let $\mathbb{X}:=\{{\mathbf{X}}^1,\dots ,{\mathbf{X}}^m\}$ be a collection of q-dimensional outputs from the m simulators that form the MME; ${\mathbf{X}}^{*}$, the true unobserved process that the simulators aim to model; and ${\mathbf{Z}}_{\mathbf{X}}$, the noisy and incomplete observation of ${\mathbf{X}}^{*}$. The model requires only two assumptions: first, that the ${\mathbf{X}}^i$ are second order exchangeable and, second, that the ${\mathbf{X}}^i$ are coexchangeable with the truth ${\mathbf{X}}^{*}$, implying $\text{cov}[{\mathbf{X}}^{*},{\mathbf{X}}^i]=\Sigma ,\forall i$. As in (3), exchangeability within $\mathbb{X}$ implies (4) $${\mathbf{X}}^i={\mathcal{M}}_{\mathbf{X}}+{\mathcal{R}}_{\mathbf{X}}^i,i=1,\dots ,m,$$(4) where ${\mathcal{M}}_{\mathbf{X}}$ is a shared mean term, the ${\mathcal{R}}_{\mathbf{X}}^i$ are the zero-mean, uncorrelated residuals of each simulator with common variance, and the ${\mathcal{M}}_{\mathbf{X}}$ and ${\mathcal{R}}_{\mathbf{X}}^i$ are uncorrelated. Coexchangeability between ${\mathbf{X}}^{*}$ and $\mathbb{X}$ implies sufficiency of ${\mathcal{M}}_{\mathbf{X}}$ for ${\mathbf{X}}^{*}$. This allows us to write (5) $${\mathbf{X}}^{*}=\mathbf{A}{\mathcal{M}}_{\mathbf{X}}+{\mathbf{U}}_{\mathbf{X}}$$(5) where A is a known matrix, and ${\mathbf{U}}_{\mathbf{X}}$ represents the ensemble discrepancy that is uncorrelated with ${\mathcal{M}}_{\mathbf{X}}$ and the ${\mathcal{R}}_{\mathbf{X}}^i$. The data, ${\mathbf{Z}}_{\mathbf{X}}$ are modeled as (6) $${\mathbf{Z}}_{\mathbf{X}}={\mathbf{H}}_{\mathbf{X}}{\mathbf{X}}^{*}+{\mathbf{W}}_{\mathbf{X}}$$(6) for measurement error ${\mathbf{W}}_{\mathbf{X}}$, and known incidence matrix ${\mathbf{H}}_{\mathbf{X}}$. The statistical model defined by (4)–(6) will hereafter by referred to as the coexchangeable model; a graphical representation is provided in Figure 1(a). Inference for this model makes use of Bayes linear sufficiency of the ensemble mean, $\overline{\mathbf{X}}=\frac{1}{m}{\sum }_{i=1}^m{\mathbf{X}}^i$, for updating by $\mathbb{X}$ and is therefore very efficient. Updated beliefs on ${\mathbf{X}}^{*}$ is done in two stages: first the update of our beliefs by the ensemble, and second by the data; the inferential procedure is outlined in the supplementary material. Climate post-processing routines, as in Steiger et al. (2014), are subsumed by the coexchangeable model, and simply describe (6) with beliefs $\mathbb{E}[{\mathbf{X}}^{*}]=\overline{\mathbf{X}}$ and $\text{var}[{\mathbf{X}}^{*}]=\mathbf{P}$. Calculations that rely on ensemble methods approximate P by the empirical covariance matrix of the ensemble, or some linear transformation thereof (see definitions in Whitaker and Hamill 2002). Our contributions, that build on the coexchangeable model, may also be considered as similar developments to such post-processing routines popular in the climate sciences.

Fig. 1 Graphical representation of the coexchangeable and coexchangeable process models. Boxes represent observed quantities, dashed circles represent unobserved quantities over which we make prior belief specifications, and solid circles represent unobserved quantities for which we calculate updated beliefs. Analogous to conditional independence in probabilistic models, arrows may be used to identify Bayes linear sufficiency between quantities. For simplicity, residual terms are omitted.

4 The Coexchangeable Process Model

4.1 Exchangeable Processes

Consider each simulator as producing output pairs of q-dimensional fields $({\mathbf{X}}^i,{\mathbf{Y}}^i)$ with the corresponding fields in reality denoted $({\mathbf{X}}^{*},{\mathbf{Y}}^{*})$, for which we have partial observations ${\mathbf{Z}}_{\mathbf{X}}$ and ${\mathbf{Z}}_{\mathbf{Y}}$, made with some error. In what follows we define SST by X and SIC by Y. Changes in SST and SIC are driven by complex physical relationships. We capture structural dependencies between random quantities $(\mathbf{X},\mathbf{Y})$, by first only imposing the coexchangeable model for ${\mathbf{X}}^i,{\mathbf{X}}^{*}$, and then considering exchangeability judgments over the processes ${\mathbf{Y}}^i$, given ${\mathbf{X}}^i$, and coexchangeability of ${\mathbf{Y}}^{*}$, given ${\mathbf{X}}^{*}$. Whilst it is enticing to consider the coexchangeable model over both fields simultaneously, it is deficient here in two ways. First, due to the dependence of ${\mathbf{Y}}^i$ on ${\mathbf{X}}^i$, the assumption of second order exchangeability between the ${\mathbf{Y}}^i$ is violated. Second, the natural way to construct the model is to define the process of ${\mathbf{Y}}^i$ given ${\mathbf{X}}^i$ parameterized by some ${\mathit{\beta }}^i$ as in a hierarchical model. Inference on the ${\mathit{\beta }}^i$ is of scientific interest; however, the coexchangeable model does not allow for this. The model sophistication required to handle the hierarchical structure requires nuanced judgments of conditional second order exchangeability. In this section, we simply state the exchangeable process model, and defer detailed discussion of the conditional exchangeability judgments to Section 4.2.

For any single simulator we model $\mathbb{E}[{\mathbf{Y}}^i]=\mathcal{M}({\mathbf{X}}^i;{\mathit{\beta }}^i)$ as a process of ${\mathbf{X}}^i$ parameterized by some ${\mathit{\beta }}^i$, specific to the ith simulator. The mean function $\mathcal{M}({\mathbf{X}}^i;{\mathit{\beta }}^i)$ may represent any relationship between ${\mathbf{X}}^i$ and ${\mathit{\beta }}^i$ and, to infer ${\mathit{\beta }}^i$ from ${\mathbf{Y}}^i$, we need only specify a joint inner product space in which they reside. We set $\mathcal{M}({\mathbf{X}}^i;{\mathit{\beta }}^i)=\varphi ({\mathbf{X}}^i){\mathit{\beta }}^i$ where $\varphi (·)$ maps ${\mathbf{X}}^i$ to some specified family of basis functions. In our application, we use a monotonic decreasing spline basis for $\varphi (·)$ to reflect the property that warmer SSTs will tend to lead to less SICs. We may then write (7) $${\mathbf{Y}}^i=\varphi ({\mathbf{X}}^i){\mathit{\beta }}^i+{\mathcal{R}}_{\mathbf{Y}}^i$$(7) where ${\mathcal{R}}_{\mathbf{Y}}^i$ is some associated residual term. A simple example to consider is where ${\mathbf{X}}^i$ and ${\mathbf{Y}}^i$ represent temporal observations at a single location, $\varphi ({\mathbf{X}}^i)={\mathbf{X}}^i$ and ${\mathit{\beta }}^i$ is a scalar. This simply models ${\mathbf{X}}^i$ and ${\mathbf{Y}}^i$ as scalar multiples of each other, where the multiplier is simulation dependent. Such examples are common when modeling emergent constraints in climate modeling. Note, the choice to model $\mathcal{M}({\mathbf{X}}^i,{\mathit{\beta }}^i)$ linear in ${\mathit{\beta }}^i$ is not required by the theory but it allows for natural specification of the inner product space and leads to sufficiency arguments, discussed below, that aid computation. To complete the analogy with a Bayesian hierarchical model, we impose second order exchangeability over ${\mathit{\beta }}^1,{\mathit{\beta }}^2,\dots$ so that (8) $${\mathit{\beta }}^i={\mathcal{M}}_{\mathit{\beta }}+{\mathcal{R}}_{\mathit{\beta }}^i,$$(8) with expectation $\mathbb{E}[{\mathit{\beta }}^i]={\mathcal{M}}_{\mathit{\beta }}$ and covariance $\text{cov}[{\mathit{\beta }}^i,{\mathit{\beta }}^{i\prime }]=\text{var}[{\mathcal{M}}_{\mathit{\beta }}]+{\mathbf{1}}_{\{i=i\prime \}}\text{var}[{\mathcal{R}}_{\mathit{\beta }}^i]$.

Conditional exchangeability does not lead to sufficiency of the sample mean and so we must calculate the belief updates in the much larger joint inner product space. Define ${\Phi }_i:=\varphi ({\mathbf{X}}^i)$; equations (7)(7) $${\mathbf{Y}}^i=\varphi ({\mathbf{X}}^i){\mathit{\beta }}^i+{\mathcal{R}}_{\mathbf{Y}}^i$$(7) and (8)(8) $${\mathit{\beta }}^i={\mathcal{M}}_{\mathit{\beta }}+{\mathcal{R}}_{\mathit{\beta }}^i,$$(8) imply that the joint space of the ${\mathbf{Y}}^i$ and the ${\mathit{\beta }}^i$ can be formed as (9) $$\left[\begin{array}{c}{\mathbf{Y}}^1\\ ⋮\\ {\mathbf{Y}}^m\\ {\mathit{\beta }}^1\\ ⋮\\ {\mathit{\beta }}^m\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}{\Phi }_1& \cdots & {\mathbf{0}}_{q\times k}\\ ⋮& \ddots & ⋮\\ {\mathbf{0}}_{q\times k}& \cdots & {\Phi }_m\end{array}& {\mathbf{0}}_{\mathit{\text{qm}}\times k}\\ {\mathbf{0}}_{\mathit{\text{km}}\times \mathit{\text{km}}}& \begin{array}{c}{\mathbf{J}}_{m\times 1}\otimes {\mathbf{I}}_k\end{array}\end{array}\right]\left[\begin{array}{c}{\mathit{\beta }}^1\\ ⋮\\ {\mathit{\beta }}^m\\ {\mathcal{M}}_{\mathit{\beta }}\end{array}\right]+\left[\begin{array}{c}{\mathcal{R}}_{\mathbf{Y}}^1\\ ⋮\\ {\mathcal{R}}_{\mathbf{Y}}^m\\ {\mathcal{R}}_{\mathit{\beta }}^1\\ ⋮\\ {\mathcal{R}}_{\mathit{\beta }}^m\end{array}\right],$$(9) where ${\mathbf{0}}_{a\times b}$ is a a × b zero matrix, ${\mathbf{J}}_{a\times b}$ is a a × b matrix of ones, and ${\mathbf{I}}_a$ is the a × a identity matrix.

Calculating the belief updates using (9) is not immediately obvious as the random quantities ${\mathit{\beta }}^i$ appear on both the left-hand side as data and the right-hand side as unknown parameters. Noting that (8) is equivalently restated as $0={\mathcal{M}}_{\mathit{\beta }}+{\mathcal{R}}_{\mathit{\beta }}^i-{\mathit{\beta }}^i$, and using the re-parameterization of Hodges (1998), we can re-express (9) as (10) $$\left[\begin{array}{c}{\mathbf{Y}}^1\\ ⋮\\ {\mathbf{Y}}^m\\ {\mathbf{0}}_{\mathit{\text{km}}\times 1}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}{\Phi }_1& \cdots & {\mathbf{0}}_{q\times k}\\ ⋮& \ddots & ⋮\\ {\mathbf{0}}_{q\times k}& \cdots & {\Phi }_m\end{array}& {\mathbf{0}}_{\mathit{\text{qm}}\times k}\\ -{\mathbf{I}}_{\mathit{\text{km}}}& \begin{array}{c}{\mathbf{J}}_{m\times 1}\otimes {\mathbf{I}}_k\end{array}\end{array}\right]\left[\begin{array}{c}{\mathit{\beta }}^1\\ ⋮\\ {\mathit{\beta }}^m\\ {\mathcal{M}}_{\mathit{\beta }}\end{array}\right]+\left[\begin{array}{c}{\mathcal{R}}_{\mathbf{Y}}^1\\ ⋮\\ {\mathcal{R}}_{\mathbf{Y}}^m\\ {\mathcal{R}}_{\mathit{\beta }}^1\\ ⋮\\ {\mathcal{R}}_{\mathit{\beta }}^m\end{array}\right]$$(10) with the familiar linear canonical form $\mathbf{Y}=\mathbf{\text{XB}}+\mathbf{E}$. As we consider the joint representation, expectation is taken jointly over the ${\mathbf{Y}}^i$ and ${\mathit{\beta }}^i$. The adjusted expectation, ${\mathbb{E}}_{\mathbf{Y}}[\mathbf{B}]$, and adjusted variance, ${\text{var}}_{\mathbf{Y}}[\mathbf{B}]$, are calculated via (1) and (2); the joint belief specifications and specific updating equations are provided in the supplementary material. We note that the specification of (9) differs from the standard approach to modeling exchangeable processes in Bayes linear statistics where inference only concerns ${\mathcal{M}}_{\mathit{\beta }}$ by, in effect, substituting (8) into (7) (see Goldstein and Wooff 1998). By jointly modeling the ${\mathit{\beta }}^i$ and ${\mathcal{M}}_{\mathit{\beta }}$ we provide a closer analogy to Bayesian hierarchical models.

As we are required to update our beliefs in the joint inner product space, computation of ${\mathbb{E}}_{\mathbf{Y}}[\mathbf{B}]$ and ${\text{var}}_{\mathbf{Y}}[\mathbf{B}]$ can be difficult. These calculations require an order $\mathcal{O}(m^3{(q+k)}^3)$ matrix inversion of $\text{var}[\mathbf{Y}]$, where m is the number of simulations in the MME, q is the dimension of the climate simulation, and k is the dimension of the ${\mathit{\beta }}^i$. Theorem 4.1 shows that a Bayes linear sufficiency argument can be made that permits a smaller computational order of $\mathcal{O}(8k^3{m}^3)$ for k < q, thus enabling efficient inference. Theorem 4.1 says, in effect, that the beliefs on B may be equivalently updated by the projection of the ${\mathbf{Y}}^i$ onto the k-dimensional column space of ${\Phi }_i$. For climate models, q is generally very large. For many basis designs $k\ll q$, and when groups, i, index MIP simulations, m is generally small.

Theorem 4.1.

Let $\widehat{\mathit{\beta }}=({\widehat{\mathit{\beta }}}^1,\dots ,{\widehat{\mathit{\beta }}}^m)$ with ${\widehat{\mathit{\beta }}}^i={({\Phi }_i^{⊺}{\Phi }_i)}^{†}{\Phi }_i^{⊺}{\mathbf{Y}}^i$. Then $\widehat{\mathit{\beta }}$ is Bayes linear sufficient for Y for adjusting B if the column space of projection matrix ${\mathrm{P}}_i={\Phi }_i{({\Phi }_i^{⊺}{\Phi }_i)}^{†}{\Phi }_i^{⊺},C({\mathrm{P}}_i)$, is invariant over i, that is, $C({\mathrm{P}}_i)=\mathcal{C},\forall i$.

Proof.

Proof available in the supplementary material. □

Matrix ${\Phi }_i$ being of full rank is sufficient, but not necessary, in satisfying the condition $C({\mathrm{P}}_i)=\mathcal{C},\forall i$ in Theorem 4.1, and most sensible basis designs will ensure this. Assuming $C({\mathrm{P}}_i)=\mathcal{C},\forall i$, and appealing to the sufficiency of $\widehat{\mathit{\beta }}$ for Y for adjusting B, we write the joint specification of the exchangeability judgments made in (10) as (11) $$\left[\begin{array}{c}{\widehat{\mathit{\beta }}}^1\\ ⋮\\ {\widehat{\mathit{\beta }}}^m\\ {\mathbf{0}}_{\mathit{\text{km}}\times 1}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{c}{\mathbf{I}}_{\mathit{\text{km}}}\end{array}& {\mathbf{0}}_{\mathit{\text{km}}\times k}\\ -{\mathbf{I}}_{\mathit{\text{km}}}& \begin{array}{c}{\mathbf{J}}_{m\times 1}\otimes {\mathbf{I}}_k\end{array}\end{array}\right]\left[\begin{array}{c}{\mathit{\beta }}^1\\ ⋮\\ {\mathit{\beta }}^m\\ {\mathcal{M}}_{\mathit{\beta }}\end{array}\right]+\left[\begin{array}{c}{\mathcal{R}}_{\widehat{\mathit{\beta }}}^1\\ ⋮\\ {\mathcal{R}}_{\widehat{\mathit{\beta }}}^m\\ {\mathcal{R}}_{\mathit{\beta }}^1\\ ⋮\\ {\mathcal{R}}_{\mathit{\beta }}^m\end{array}\right]$$(11) where ${\mathcal{R}}_{\widehat{\mathit{\beta }}}^i={({\Phi }_i^{⊺}{\Phi }_i)}^{†}{\Phi }_i^{⊺}{\mathcal{R}}_{\mathbf{Y}}^i$. Calculation of ${\mathbb{E}}_{\widehat{\mathit{\beta }}}[\mathbf{B}]$ and ${\text{var}}_{\widehat{\mathit{\beta }}}[\mathbf{B}]$ follows (1) and (2); the specific equations are given in the supplementary material. We also provide partial updating equations for ${\mathbb{E}}_{\widehat{\mathit{\beta }}}[{\mathcal{M}}_{\mathit{\beta }}]$ and ${\text{var}}_{\widehat{\mathit{\beta }}}[{\mathcal{M}}_{\mathit{\beta }}]$ which are calculated with computational order $\mathcal{O}(k^3{m}^3)$ and may be used in equations (S4) and (S5) that update the full coexchangeable process model.

4.2 Repeated Observations of a Process

Assume, within each simulator output i, we observe a sequence $({\mathbf{X}}_1^i,{\mathbf{Y}}_1^i),({\mathbf{X}}_2^i,{\mathbf{Y}}_2^i),\dots$ of observation/covariate pairings; here, each $({\mathbf{X}}_t^i,{\mathbf{Y}}_t^i)$ is some p-dimensional sub-level process contained in $({\mathbf{X}}^i,{\mathbf{Y}}^i)$ and the ${\mathit{\beta }}^i$ are invariant to the dimension indexed by t. In the simple example provided above, each $({\mathbf{X}}_t^i,{\mathbf{Y}}_t^i)$ are pairings of scalar observations, and so p = 1, and ${\mathit{\beta }}^i$ is invariant in time. Similarly, for our application, each t indexes time, but the ${\mathbf{X}}_t^i$ and ${\mathbf{Y}}_t^i$ are spatial observations within the spatio-temporal ${\mathbf{X}}^i$ and ${\mathbf{Y}}^i$; though we could instead index space, both space and time, or some other feature of the process. The exchangeable process model implicitly requires an assumption of conditional exchangeability in (7), the second order equivalent of the exchangeability judgments used in Bayesian regression (for discussion see Williamson and Sansom 2019). Specifically, conditional second order exchangeability implies that the ${\mathbf{Y}}_1^i,{\mathbf{Y}}_2^i,\dots$ are second order exchangeable given some ${\mathbf{X}}_t^i$ that is constant for all t. We may thus specify a mean process, $\mathbb{E}[{\mathbf{Y}}_t^i]=\mathcal{M}({\mathbf{X}}_t^i;{\mathit{\beta }}^i)$ for some known covariate ${\mathbf{X}}_t^i$ and unknown parameter ${\mathit{\beta }}^i$; above, we assume the mean process to be $\mathbb{E}[{\mathbf{Y}}_t^i]=\varphi ({\mathbf{X}}_t^i){\mathit{\beta }}^i$. According to the second order representation theorem, (12) $${\mathbf{Y}}_t^i=\varphi ({\mathbf{X}}_t^i){\mathit{\beta }}^i+{\mathcal{R}}_{{\mathbf{Y}}_t}^i$$(12) where $\text{cov}[{\mathbf{Y}}_t^i,{\mathbf{Y}}_{t\prime }^i]=\text{cov}[\mathcal{M}({\mathbf{X}}_t^i;{\mathit{\beta }}^i),\mathcal{M}({\mathbf{X}}_{t\prime }^i;{\mathit{\beta }}^i)]+{\mathbf{1}}_{\{t=t\prime \}}$ $\text{var}[{\mathcal{R}}_{{\mathbf{Y}}_t}^i]$ and we arrive at (7) by building the joint representation, that is, stacking the instances of t. In can be intuitive to think of the full simulator outputs $({\mathbf{X}}^i,{\mathbf{Y}}^i)$ as matrices with the variant and invariant dimensions (here, space and time) indexed across the rows and columns, respectively. The mathematics in Section 4.1 simply then require ${\mathbf{X}}^i$ and ${\mathbf{Y}}^i$ to be substituted by their vectorized equivalents, $\text{vec}({\mathbf{X}}^i)$ and $\text{vec}({\mathbf{Y}}^i)$.

4.3 Reality as a Coexchangeable Process

We now show how judgments of coexchangeability may be made to incorporate (10) or (11) into the coexchangeable model. Define the true process of ${\mathbf{Y}}^i$ that the simulators attempt to resolve as ${\mathbf{Y}}^{*}$. An assumption coexchangeability of ${\mathbf{Y}}^{*}$ and ${\mathbf{Y}}^i$ given fixed ${\mathbf{X}}^i=\mathbf{X}$ $\forall i$, and hence ${\Phi }_i=\Phi$ $\forall i$, is equivalent to assuming coexchangeability of ${\mathbf{Y}}^{*}$ and the ${\mathit{\beta }}^i$. Bayes linear sufficiency of ${\mathcal{M}}_{\mathit{\beta }}$ for the ${\mathbf{Y}}^i$ for adjusting ${\mathbf{Y}}^{*}$ follows. We write a model for ${\mathbf{Y}}^{*}$ such that (13) $${\mathbf{Y}}^{*}={\mathbf{A}}_{\mathbf{Y}}{\mathcal{M}}_{\mathit{\beta }}+{\mathbf{U}}_{\mathbf{Y}},$$(13) where ${\mathbf{U}}_{\mathbf{Y}}$ is a model-mismatch term uncorrelated with the ${\mathbf{Y}}^i$ and the assumption of coexchangeability permits ${\mathbf{A}}_{\mathbf{Y}}$ to be any matrix of suitable dimensions. The choice of ${\mathbf{A}}_{\mathbf{Y}}=\varphi ({\mathbf{X}}^{*})$ is obvious in our context but different choices are permissible should they make sense to the application. Finally, the data, ${\mathbf{Z}}_{\mathbf{Y}}$ are modeled as (14) $${\mathbf{Z}}_{\mathbf{Y}}={\mathbf{H}}_{\mathbf{Y}}{\mathbf{Y}}^{*}+{\mathbf{W}}_{\mathbf{Y}},$$(14) for measurement error ${\mathbf{W}}_{\mathbf{Y}}$, and known incidence matrix ${\mathbf{H}}_{\mathbf{Y}}$.

Similar to Rougier, Goldstein, and House (2013), we assume sufficiency of ${\mathbf{Y}}^{*}$ for ${\mathbf{Z}}_{\mathbf{Y}}$, allowing us to update beliefs on ${\mathbf{Y}}^{*}$ in two stages. The assumption of sufficiency of ${\mathbf{Y}}^{*}$ for ${\mathbf{Z}}_{\mathbf{Y}}$ is akin to an assumption of conditional independence between ${\mathbf{Z}}_{\mathbf{Y}}$ and the ${\mathbf{Y}}^i$ in a probabilistic analysis. The two-stage update of beliefs on ${\mathbf{Y}}^{*}$ may thus be equivalently thought of as either the joint update by the ${\mathbf{Y}}^i$ and ${\mathbf{Z}}_{\mathbf{Y}}$ or as sequential updates by the ${\mathbf{Y}}^i$ and then ${\mathbf{Z}}_{\mathbf{Y}}$. As with the coexchangeable model in Section 3.3, we show an graphical representation of coexchangeable process model in Figure 1(b); the equations for the two-stage belief update are provided in the supplementary material.

5 Joint Reconstruction of Paleo Sea-Surface Temperature and Sea-Ice with PMIP3 and PMIP4

5.1 Paleodata

Geological, ecological and geochemical measurements of the LGM have large associated uncertainties, and these uncertainties are further compounded by relating the measured processes into proxies for SST and SIC. We use judgments from subject matter experts to account for sampling bias and the errors in the proxy-data, which are considered to be systematic in space. In some cases we directly incorporate these judgments into the belief structure of the model, namely, expectations, variances, and covariances. In other cases, similar to Rougier et al. (2022), we use “pseudo-observations” to reflect subject matter experts’ judgments in sparsely observed areas.

SST data, ${\mathbf{Z}}_{\mathbf{X}}$, uses proxy measurements obtained from sea-bed sediment core samples recorded either as annual or summer means. The measurements are shown in Figure 2; annual and summer means are depicted with points and triangles, respectively. There is very likely some strong observation bias in the foraminifera-based proxy measurements from the Arctic and Nordic seas, approximately north of ${62}^{°}$ N. In this region, very cold water and full sea-ice coverage are common; both inhospitable conditions for most foraminifera species. The nonexistence of foraminifera in ocean sediments are rarely reported, since the absence of foraminifera leaves little to be analyzed in ecological or geochemical studies. Thus, observations may be biased toward the warm climate events naturally present within the inter-decadal variability, when foraminifera are found. We account for this through the specification of the measurement error term in (6), ${\mathbf{W}}_{\mathbf{X}}$. Define ${\mathrm{B}}_{\mathrm{N}}$ as a set of indices that index observations from the Nordic Seas in ${\mathbf{Z}}_{\mathbf{X}}$. We set $\mathbb{E}{[{\mathbf{W}}_{\mathbf{X}}]}_b=2$ and $\mathbb{E}{[{\mathbf{W}}_{\mathbf{X}}]}_{b^{*}}=0$ for $b\in {\mathrm{B}}_{\mathrm{N}}$ and $b^{*}\notin {\mathrm{B}}_{\mathrm{N}}$. As the bias originates from a systematic reporting error we believe the measurement errors to be correlated. We partition $\text{var}[{\mathbf{W}}_{\mathbf{X}}]={\mathrm{V}}_{\mathbf{D}}(\mathbf{I}+{\mathrm{V}}_{\mathbf{B}}){\mathrm{V}}_{\mathbf{D}}$ and set ${\mathrm{V}}_{\mathbf{D}}$ as the diagonal matrix of the reported standard deviations in the MARGO dataset, and ${\mathrm{V}}_{\mathbf{B}[b,b\prime ]}=1\forall b,b\prime \in {\mathrm{B}}_{\mathrm{N}}$ where $b\ne b\prime$, and 0 otherwise. This represents the weakest possible belief specification as it makes the measurement error in the Nordic seas perfectly correlated, in essence reducing the information into a single observation. Without accounting for the bias in these observations, the Nordic Seas would be too warm to support sea-ice, which is known from geological records not to be the case. The specification of $\mathbb{E}{[{\mathbf{W}}_{\mathbf{X}}]}_b=2$ comes from subject matter experts; even still, due to the imposed correlation structure the analysis is not sensitive to this specification. Finally, we set incidence matrix in (6), ${\mathbf{H}}_{\mathbf{X}}={\mathbf{H}}_{\mathbf{X}}^{\mathrm{T}}\otimes {\mathbf{H}}_{\mathbf{X}}^{\mathrm{S}}$, where ${\mathbf{H}}_{\mathbf{X}}^{\mathrm{T}}$ calculates either the annual or summer means of ${\mathbf{X}}^{*}$ as necessary, and ${\mathbf{H}}_{\mathbf{X}}^{\mathrm{S}}$ spatially interpolates these averages from simulator’s spatial grid to the data locations.

Fig. 2 Proxy-based measurements of LGM SST and maximum sea-ice extents. Annual and summer mean SST’s are represented by points and triangles, respectively. Sea-ice extents are represented by a red dashed line. The Southern sea-ice extent is as reported in Gersonde et al. (2005) and the Northern sea-ice extent is provided by coauthors.

Point-wise proxy measurements of SIC are difficult to interpret and unreliable. More robust are estimates of maximum sea-ice extent. The Northern and Southern extents used for ${\mathbf{Z}}_{\mathbf{Y}}$ are shown in Figure 2 by the red dashed lines. ${\mathbf{Z}}_{\mathbf{Y}}$ is a 4145-dimensional vector, that is, the same spatial resolution as the numerical models, that records a 1 at spatial locations within the extent boundaries, and a 0 outside. Incidence matrix, ${\mathbf{H}}_{\mathbf{Y}}$ in (14) is a p × q matrix that collates, from ${\mathbf{Y}}^{*}$, the February SIC from the Northern hemisphere and August SIC from the Southern hemisphere. We specify measurement error, ${\mathbf{W}}_{\mathbf{Y}}$, to be spatially correlated and certain of SIC estimates in the poles, where we are confident that there is full sea-ice coverage, and mid-latitude and equatorial regions, where we are confident there is zero sea-ice. We set $\text{var}[{\mathbf{W}}_{\mathbf{Y}}]=K\text{cor}[{\mathbf{W}}_{\mathbf{Y}}]K$ where K is a diagonal matrix that represents our marginal uncertainty and $\text{cor}[{\mathbf{W}}_{\mathbf{Y}}]$ is spatially correlated error; descriptions of K and $\text{cor}[{\mathbf{W}}_{\mathbf{Y}}]$ are given in the supplementary material.

5.2 Fitting the Coexchangeable Process Model

To fit the coexchangeable process model we first model SST via the coexchangeable model described in Section 3.3 and the process of SIC given SST is modeled using the methodology developed in Section 4. As described above, we build a MME by selecting a single representative simulation from each of the PMIP3 and PMIP4 modeling groups, with the exception of the HadCM3 model simulations where we make use of all three available PMIP4 simulations that use different ice sheet boundary conditions. The m = 13 models selected are (15) $$\begin{array}{cc}\mathcal{S}=\{\hfill & \mathtt{\text{CNRM}}-\mathtt{3},\mathtt{\text{IPSL}}-\mathtt{3},\mathtt{\text{MIROC}}-\mathtt{3},\mathtt{\text{MPI}}-\mathtt{3},\hfill \\ \hfill & \mathtt{\text{CCSM}}\mathtt{4}-\mathtt{3},\mathtt{\text{GISS}}-\mathtt{3},\mathtt{\text{MRI}}-\mathtt{3},\mathtt{\text{AWI}}-\mathtt{4},\hfill \\ \hfill & \mathtt{\text{MIROC}}-\mathtt{4},\mathtt{\text{MPI}}-\mathtt{4},\mathtt{\text{HadCM}}\mathtt{3}-\mathtt{\text{PMIP}},\hfill \\ \hfill & \mathtt{\text{HadCM}}\mathtt{3}-\mathtt{\text{Glac}}\mathtt{1D},\mathtt{\text{HadCM}}\mathtt{3}-\mathtt{\text{Ice}}\mathtt{6G}\}.\hfill \end{array}$$(15)

The assumption of coexchangeability is a prior judgment, and at the time of the analysis each of these models was deemed coexchangeable by the project’s SMEs. Each ensemble member was projected onto the FAMOUS ocean grid with 4145 spatial locations, and we use the 12 monthly means of SST and SIC.

5.2.1 Sea-Surface Temperature

Following notation in Section 3.3, define the ensemble SSTs as $\mathbb{X}=\{{\mathbf{X}}^1,\dots ,{\mathbf{X}}^m\}$, and assume second order exchangeability within the ensemble. This leads to the representation in (4), for which we require prior specification of $\text{var}[{\mathcal{M}}_{\mathbf{X}}]$ and $\text{var}[{\mathcal{R}}_{\mathbf{X}}^i]$. Coexchangeabilty of ${\mathbf{X}}^i$ and ${\mathbf{X}}^{*}$ leads to (5), for which we require prior specification of model mismatch terms $\mathbb{E}[{\mathbf{U}}_{\mathbf{X}}],\text{ var}[{\mathbf{U}}_{\mathbf{X}}]$ and incidence matrix A. We specify $\text{var}[{\mathcal{R}}_{\mathbf{X}}]={\alpha }^2\text{var}[{\mathcal{M}}_{\mathbf{X}}]$ so that $\text{var}[\mathbf{X}]=(1+{\alpha }^2)\text{var}[{\mathcal{M}}_{\mathbf{X}}]$ and, following Rougier, Goldstein, and House (2013), set $\mathbb{E}[{\mathbf{U}}_{\mathbf{X}}]=\mathbf{0}$ and $\mathbf{A}=\mathbf{I}$.

For most climate fields, including SST, we can exploit spatio-temporal structure so that computation of adjusted beliefs is scalable. The most obvious way to do this is to assume separability through space and time so that $\text{var}[\mathbf{X}]$ and $\text{var}[{\mathbf{X}}^{*}]$, and hence $\text{var}[{\mathbf{U}}_{\mathbf{X}}]$, have a Kronecker structure. For example, $\text{var}[\mathbf{X}]=\text{var}[{\mathbf{X}}_{\mathrm{T}}]\otimes \text{var}[{\mathbf{X}}_{\mathrm{S}}]$ where $\mathrm{T}$ denotes time and $\mathrm{S}$ denotes space. If we similarly partition $\text{var}[{\mathbf{U}}_{\mathbf{X}}]=\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{T}}}]\otimes \text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{S}}}]$ and equate either $\text{var}[{\mathbf{X}}_{\mathrm{T}}]=\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{T}}}]$ or $\text{var}[{\mathbf{X}}_{\mathrm{S}}]=\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{S}}}]$ computation of our adjusted beliefs is efficient. Note that this assumption is weaker than the application in Rougier, Goldstein, and House (2013) where it is assumed that ${\kappa }^2\text{var}[\mathbf{X}]=\text{var}[\mathbf{U}]$. Here, we set $\text{var}[{\mathbf{X}}_{\mathrm{T}}]=\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{T}}}]={\mathbf{J}}_{n\times n}$ as our subject matter experts believe that discrepancies between the ${\mathbf{X}}^i$ and reality are predominantly spatial and constant in time.

We set $\text{var}[{\mathbf{X}}_{\mathrm{S}}]$ as the positive semidefinite matrix that minimizes the distance between $\text{var}[\mathbf{X}]$ and the empirical covariance matrix of $\mathbb{X},{\mathbf{S}}_{\mathbf{X}}$, under the Frobenius norm. This specification follows similar arguments to Rougier, Goldstein, and House (2013) where $\text{var}[\mathbf{X}]={\mathbf{S}}_{\mathbf{X}}$, but preserves Kronecker structure in $\text{var}[\mathbf{X}]$ to allow for efficient computation. Details of this calculation are in the supplementary material. Elements of $\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{S}}}]$ are defined via the C⁴-Wendland covariance function such that (16) $$\begin{array}{c}\text{var}{\left[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{S}}}\right]}_{[s,s\prime ]}={\kappa }^2\left(1+\tau \frac{d(s,s\prime )}{c}+\frac{{\tau }^2-1}{3}\frac{d{(s,s\prime )}^2}{c^2}\right)\\ {\left(1-\frac{d(s,s\prime )}{c}\right)}_{+}^{\tau }\end{array}$$(16) where $\tau \ge 6,c\in (0,\pi ],{(a)}_{+}=\mathrm{\text{max}}(0,a)$, and d(i, j) is the geodesic distance between locations i and j. The C⁴-Wendland covariance function is commonly chosen so as to define a smooth process on the sphere; (see, e.g., Astfalck et al. 2019). Parameters are specified as $\kappa =1.61$, c = 0.92, and τ = 6; these values are selected to represent the subject matter experts’ beliefs as to the magnitude and correlation lengths of ${\mathbf{U}}_{\mathbf{X}}$. A sensitivity analysis is provided in the supplementary material to highlight how these judgments influence inference.

The two-stage Bayes linear update follows Rougier, Goldstein, and House (2013). As with Rougier, Goldstein, and House (2013), we assume the first update ${\mathbb{E}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]$ is well approximated by ${\mathbb{E}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]\approx \overline{\mathbf{X}}$ and we calculate ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]=\left(\frac{{\alpha }^2}{m+{\alpha }^2}\right)\text{var}[{\mathcal{M}}_{\mathbf{X}}]+\text{var}[{\mathbf{U}}_{\mathbf{X}}]$ so ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]\to \text{var}[\mathbf{U}]$ as $m\to \infty$. Here, we choose ${\alpha }^2=1$ and have m = 13 and so do not assume ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]\approx \text{var}[\mathbf{U}]$. From specifications $\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{T}}}]=\text{var}[{\mathbf{X}}_{\mathrm{T}}]$ and $\text{var}[\mathbf{X}]=2\text{var}[{\mathcal{M}}_{\mathbf{X}}],{\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]$ has Kronecker structure so that ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]={\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}_{\mathrm{T}}^{*}]\otimes {\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}_{\mathrm{S}}^{*}]$, where ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}_{\mathrm{T}}^{*}]={\mathbf{J}}_{n\times n}$ and ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}_{\mathrm{S}}^{*}]=\frac{\text{var}[{\mathbf{X}}_{\mathrm{S}}]}{2(m+1)}+\text{var}[{\mathbf{U}}_{{\mathbf{X}}_{\mathrm{S}}}]$. The above specifications lead to a prior predictive for ${\mathbf{X}}^{*}$ that is warmer than our true beliefs at the poles where we are certain there was full sea ice coverage and so the SST must be -1.92 °C. This is problematic due to the data sparsity in the poles, and so we correct our prior by adding 10 equally longitudinally spaced pseudo-observations at 80°N and 80°S, each of -1.92 °C. Figure 3(a) and (b) plot ${\mathbb{E}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]$ and marginal standard deviation of ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]$, respectively, for January.

Fig. 3 Adjusted beliefs of January SST: (a) expectation of SST adjusted by $\overline{\mathbf{X}}$, equal to the ensemble mean; (b) marginal standard deviation of SST adjusted by $\overline{\mathbf{X}}$; (c) expectation of SST adjusted by $\overline{\mathbf{X}}$ and ${\mathbf{Z}}_{\mathbf{X}}$; (d) marginal standard deviation of SST adjusted by $\overline{\mathbf{X}}$ and ${\mathbf{Z}}_{\mathbf{X}}$; (e) contribution of the data to our expected beliefs of SST ${\mathbb{E}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}[{\mathbf{X}}^{*}]-\overline{\mathbf{X}}$; and (f), marginal standard deviation of the ensemble. All plots are shown in degrees Celsius.

The second update is calculated by (17) $$\begin{array}{c}{\mathbb{E}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}\left[{\mathbf{X}}^{*}\right]=\overline{\mathbf{X}}+{\text{var}}_{\overline{\mathbf{X}}}\left[{\mathbf{X}}^{*}\right]{\mathbf{H}}_{\mathbf{X}}^{⊺}{\text{var}}_{\overline{\mathbf{X}}}{\left[{\mathbf{Z}}_{\mathbf{X}}\right]}^{-1}\\ ({\mathbf{Z}}_{\mathbf{X}}-{\mathbf{H}}_{\mathbf{X}}\overline{\mathbf{X}}-\mathbb{E}[{\mathbf{W}}_{\mathbf{X}}])\text{, and}\end{array}$$(17) (18) $$\begin{array}{c}{\text{var}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}\left[{\mathbf{X}}^{*}\right]={\mathbf{J}}_{n\times n}\otimes ({\text{var}}_{\overline{\mathbf{X}}}\left[{\mathbf{X}}_{\mathrm{S}}^{*}\right]-{\text{var}}_{\overline{\mathbf{X}}}\left[{\mathbf{X}}_{\mathrm{S}}^{*}\right]{\left({\mathbf{H}}_{\mathbf{X}}^{\mathrm{S}}\right)}^{⊺}\\ {\text{var}}_{\overline{\mathbf{X}}}{\left[{\mathbf{Z}}_{\mathbf{X}}\right]}^{-1}{\mathbf{H}}_{\mathbf{X}}^{\mathrm{S}}{\text{var}}_{\overline{\mathbf{X}}}\left[{\mathbf{X}}_{\mathrm{S}}^{*}\right]),\end{array}$$(18) where ${\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{Z}}_{\mathbf{X}}]={\mathbf{H}}_{\mathbf{X}}{\text{var}}_{\overline{\mathbf{X}}}[{\mathbf{X}}^{*}]{\mathbf{H}}_{\mathbf{X}}^{⊺}+\text{var}[{\mathbf{W}}_{\mathbf{X}}]$. Figure 3(c) and (d) show these updates, also for January, respectively. Figure 3(e) and (f) give indication of what information is gained from the data: Figure 3(e) shows the difference the ensemble mean and ${\mathbb{E}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}[{\mathbf{X}}^{*}]$; Figure 3(f) shows the empirical standard deviation of the ensemble, which when compared to Figure 3(d) is shown to be more uncertain. Regions of cooling are apparent westward of large continental masses, that is on the North American Pacific Coast, or the African Atlantic Coast. This is indicative of models not capturing up-welling phenomena, and is a pattern previously observed in the error between models and pre-industrial simulations (Eyring et al. 2019), as well as the cooling in the Southern Ocean, and warming the Indian and Pacific oceans. Larger uncertainty is seen in the Pacific where measurements are sparse.

5.2.2 Sea-Ice Concentration Given Sea-Surface Temperature

Following notation in Section 4 we now consider the process of SIC given SST. Define the jointly-observed ensemble of SST and SIC as $(\mathbb{X},\mathbb{Y})=\{({\mathbf{X}}^1,{\mathbf{Y}}^1),\dots ,({\mathbf{X}}^m,{\mathbf{Y}}^m)\}$ where each ${\mathbf{Y}}^i$ is a spatio-temporal vector of SIC that we model dependent on SST, ${\mathbf{X}}^i$. The ${\mathbf{Y}}^i$ comprise 12 (i.e., monthly) 4145-dimensional conditionally second order exchangeable spatial processes, ${\mathbf{Y}}_t^i$, where i indexes the ensemble member and t indexes time. We assume $\mathbb{E}[{\mathbf{Y}}_t^i]=\varphi ({\mathbf{X}}_t^i){\mathit{\beta }}^i$ which leads to the representation in (7), for which we require specifications of $\varphi (·),\mathbb{E}[{\mathit{\beta }}^i],\text{ var}[{\mathit{\beta }}^i]$, and $\text{var}[{\mathcal{R}}_{{\mathbf{Y}}_t}^i]$. We assume the ${\mathit{\beta }}^i$ to be second order exchangeable which leads to the representation in (8) that requires specification of $\mathbb{E}[{\mathcal{M}}_{\mathit{\beta }}],\text{ var}[{\mathcal{M}}_{\mathit{\beta }}]$ and $\text{var}[{\mathcal{R}}_{\mathit{\beta }}]$. Together, these exchangeability judgments lead to the representation for the ${\mathbf{Y}}^i$ in (9), and equivalently (10).

For modeling SIC given SST we require spatial variation in the physics of the process. For example, the relationship between SIC and SST is different in the Nordic seas where sea-ice is supported at warmer SSTs than other locations. To account for this, we specify $\varphi ({\mathbf{X}}_t^i)=\Psi ({\mathbf{X}}_t^i){\mathbf{\Theta }}^i$ where $\Psi ({\mathbf{Y}}_t^i)$ models the behavior of SIC and SST at individual locations using spline bases, and ${\mathbf{\Theta }}^i$ is a fixed-rank spatial basis of the spline coefficients. Note that we can model the spline coefficients individually in space, in which case $\varphi ({\mathbf{X}}_t^i)=\Psi ({\mathbf{X}}_t^i)$, but as the size of ${\mathit{\beta }}^i$ then scales with the spatial resolution; this is not feasible for climate models. Here we specify $\Psi ({\mathbf{X}}_t^i)$ with I-spline bases, a basis family commonly used for monotone functions (Ramsay 1988). To ease computation, we project the spatial basis ${\mathbf{\Theta }}^i$ onto a principal component design calculated using a projection of the ${\mathbf{Y}}^i$ onto the column space of the $\Psi ({\mathbf{X}}^i),{\widehat{\Theta }}^i$. Approximating the spatial coefficients using principal components restricts inference for the ${\mathit{\beta }}^i$ and ${\mathcal{M}}_{\mathit{\beta }}$ to linear combinations of ${\widehat{\Theta }}^i$; we argue this is appropriate for modeling the mean MME process. We use a more flexible parametrization in the model discrepancy below so that inference on ${\mathbf{Y}}^{*}$ given the SIC data is not restricted to linear combinations of the principal component design. Full specification of $\Psi ({\mathbf{X}}_t^i)$ and the ${\mathbf{\Theta }}^i$ are in the supplementary material. All ${\Phi }_i=\varphi ({\mathbf{X}}^i)$ are full rank and so the conditions of Theorem 4.1 are met and $\widehat{\mathit{\beta }}=({\widehat{\mathit{\beta }}}^1,\dots ,{\widehat{\mathit{\beta }}}^m)$ is sufficient for $(\mathbb{X},\mathbb{Y})$ for updating $\mathbf{B}=({\mathit{\beta }}^1,\dots ,{\mathit{\beta }}^m,{\mathcal{M}}_{\mathit{\beta }})$. We specify $\text{var}[{\mathcal{R}}_{{\mathbf{Y}}_t}^i]$ as a heteroscdastic error process; smaller variance is specified for very cold SSTs, where we are confident there is sea-ice, and warm SSTs, where we are confident there is no sea-ice, and larger variance is specified for SSTs approximately between -1 °C and 3 °C where sea-ice behavior is variable. As above we assume $\mathbb{E}[{\mathcal{M}}_{\mathit{\beta }}]$ is well approximated by the empirical mean of the MME members so that $\mathbb{E}[{\mathcal{M}}_{\mathit{\beta }}]=\frac{1}{m}{\sum }_{i=1}^m{\widehat{\mathit{\beta }}}^i$; further, similar to Rougier, Goldstein, and House (2013), we set $\text{var}[{\mathit{\beta }}^i]=\lambda$ where λ is a diagonal matrix of the eigenvalues calculated from the principal component decomposition described above. As before we set $\text{var}[{\mathcal{R}}_{\mathit{\beta }}]=\text{var}[{\mathcal{M}}_{\mathit{\beta }}]$ and so $\text{var}[{\mathcal{R}}_{\mathit{\beta }}]=\text{var}[{\mathcal{M}}_{\mathit{\beta }}]=\lambda /2$.

Adjusted beliefs ${\mathbb{E}}_{\widehat{\mathit{\beta }}}[\mathbf{B}]$ and ${\text{var}}_{\widehat{\mathit{\beta }}}[\mathbf{B}]$ are calculated as in the supplementary material. As mentioned in Section 4.1 an advantage of specifying our model jointly over the ${\mathit{\beta }}^i$ and ${\mathcal{M}}_{\mathit{\beta }}$ is that we obtain inference for each ${\mathit{\beta }}^i$ (as opposed to only ${\mathcal{M}}_{\mathit{\beta }}$ as in Goldstein and Wooff 1998). We look at the individual fits in Figure 4 where we plot, at four locations, the fit ensemble members prior to the coexchangeable adjustment. Location A shows a point in the Arctic where the models largely agree. Location B shows a region that occasionally sees small concentrations of sea-ice in the models but is predominantly too warm to support full sea-ice coverage. Location C shows three MME members whose relationship of SIC given SST disagree with the rest of the ensemble. Finally, location D shows sea-ice behavior at the Antarctic sea-ice edge. Given the SSTs that are observed, the differing physics predominantly lead to differences in Winter sea-ice.

Fig. 4 Spline fits, parameterized by ${\mathit{\beta }}^i$, of SIC given SST to each ensemble member prior to the coexchangeable adjustment. Shaded regions denote ±2 standard deviations.

Coexchangeability of ${\mathbf{Y}}^{*}$ and ${\mathit{\beta }}^i$ leads to (13) for which specification of $A({\mathbf{X}}^{*}),\mathbb{E}[{\mathbf{U}}_{\mathbf{Y}}]$ and $\text{var}[{\mathbf{U}}_{\mathbf{Y}}]$ is required. We set $A({\mathbf{X}}^{*})=\varphi ({\mathbf{X}}^{*})$ and model discrepancy so that ${\mathbf{U}}_{\mathbf{Y}}=\Psi ({\mathbf{X}}^{*}){\mathbf{U}}_{\Theta }$, and thus $\mathbb{E}[{\mathbf{U}}_{\mathbf{Y}}]=\Psi ({\mathbf{X}}^{*})\mathbb{E}[{\mathbf{U}}_{\Theta }]$ and $\text{var}[{\mathbf{U}}_{\mathbf{Y}}]=\Psi ({\mathbf{X}}^{*})\text{var}[{\mathbf{U}}_{\Theta }]\Psi {({\mathbf{X}}^{*})}^{⊺}$. The matrix ${\mathbf{U}}_{\Theta }$ represents the discrepancy of the spatial spline coefficients for the reality model. Should our MME contain models with high spatial resolution we could specify ${\mathbf{U}}_{\Theta }$ as a lower-rank process (e.g., with fixed-rank methods such as Cressie and Johannesson 2008); we do not find the need to do so here. We set $\mathbb{E}[{\mathbf{U}}_{\Theta }]=0$ and $\text{var}[{\mathbf{U}}_{\Theta }]={\mathbf{I}}_{l\times l}\otimes \text{var}[{\mathbf{U}}_{{\mathbf{Y}}_{\mathrm{S}}}]$ where l is the number of spline coefficients at each location and $\text{var}[{\mathbf{U}}_{{\mathbf{Y}}_{\mathrm{S}}}]$ is calculated by the C⁴-Wendland function in (16) with parameters $\kappa =0.3$, c = 4, and τ = 6; again, we include a sensitivity analysis in the supplementary material to examine the sensitivity of our inference to differing parameterizations. The first update ${\mathbb{E}}_{\widehat{\mathit{\beta }}}[{\mathbf{Y}}^{*}]$ and ${\text{var}}_{\widehat{\mathit{\beta }}}[{\mathbf{Y}}^{*}]$ is calculated given ${\mathbf{X}}^{*}$ using ${\mathbb{E}}_{\widehat{\mathit{\beta }}}[{\mathcal{M}}_{\mathit{\beta }}]$ and ${\text{var}}_{\widehat{\mathit{\beta }}}[{\mathcal{M}}_{\mathit{\beta }}]$ as previously calculated. The second update ${\mathbb{E}}_{\widehat{\mathit{\beta }},{\mathbf{Z}}_{\mathbf{Y}}}[{\mathbf{Y}}^{*}]$ and ${\text{var}}_{\widehat{\mathit{\beta }},{\mathbf{Z}}_{\mathbf{Y}}}[{\mathbf{Y}}^{*}]$ are calculated given ${\mathbf{X}}^{*}$ as in (S6) and (S7). Expectations and marginal standard deviations of both updates, calculated with ${\mathbf{X}}^{*}={\mathbb{E}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}[{\mathbf{X}}^{*}]$, are shown in Figure 5. The expected SSTs do not produce adequate sea-ice when only considering learnt relations from the MME (Figure 5(a)) and marginal standard deviation of SIC is large in regions where the expected SST is cold enough to guarantee sea-ice coverage (Figure 5(b)). Updating SIC using the data leads to more extensive sea-ice cover (Figure 5(c)) and a reduction of standard deviation everywhere except the sea-ice boundary (Figure 5(d)). We also show, in Figure 5(e) and (f) the empirical ensemble mean and standard deviation, respectively. Our final reconstruction produces more sea-ice in the Southern Hemisphere, but more crucially, drastically reduces the sea-ice uncertainty in the sea-ice interior. Similar behavior is seen in the Northern Hemisphere winter.

Fig. 5 Adjusted beliefs of August SIC given expected SST, ${\mathbf{X}}^{*}={\mathbb{E}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}[{\mathbf{X}}^{*}]$: (a) expectation of SIC adjusted by $\widehat{\mathit{\beta }}$; (b) marginal standard deviation of SIC adjusted by $\widehat{\mathit{\beta }}$; (c) expectation of SIC adjusted by $\widehat{\mathit{\beta }}$ and ${\mathbf{Z}}_{\mathbf{Y}}$; (d) marginal standard deviation of SIC adjusted by $\widehat{\mathit{\beta }}$ and ${\mathbf{Z}}_{\mathbf{Y}}$; (e) ensemble mean; and (f), marginal standard deviation of the ensemble. Values are of sea-ice concentration measured between 0 and 1.

6 Sampled Boundary Conditions and their Influence on Glacial Ice Sheet Modeling Outputs

The remit of this work was to reconstruct, with uncertainty, joint SST and SIC fields to act as boundary conditions into FAMOUS-Ice (Smith, George, and Gregory 2021), an ‘atmosphere only’ global climate model coupled to a ice sheet model. To examine the effects that the reconstruction has on the atmosphere-ice sheet model outputs we run a small ensemble varying only the SST and SIC boundary conditions. Bayes linear analysis updates our beliefs of the first two moments of ${\mathbf{X}}^{*}$ and ${\mathbf{Y}}^{*}$. Define the Cholesky decompositions ${\mathbf{L}}_{\mathbf{X}}$ and ${\mathbf{L}}_{\mathbf{Y}}$, so that ${\mathbf{L}}_{\mathbf{X}}{\mathbf{L}}_{\mathbf{X}}^{⊺}={\text{var}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}[{\mathbf{X}}^{*}]$ and ${\mathbf{L}}_{\mathbf{Y}}{\mathbf{L}}_{\mathbf{Y}}^{⊺}={\text{var}}_{\widehat{\mathit{\beta }},{\mathbf{Z}}_{\mathbf{Y}}}[{\mathbf{Y}}^{*}]$. We may probabilistically sample (19) $${\tilde{\mathbf{X}}}^{*}\sim {\mathbb{E}}_{\overline{\mathbf{X}},{\mathbf{Z}}_{\mathbf{X}}}\left[{\mathbf{X}}^{*}\right]+{\mathbf{L}}_{\mathbf{X}}\mathbf{Z}$$(19) (20) $${\tilde{\mathbf{Y}}}^{*}\sim {\mathbb{E}}_{\widehat{\mathit{\beta }},{\mathbf{Z}}_{\mathbf{Y}}}\left[{\mathbf{Y}}^{*}\right]+{\mathbf{L}}_{\mathbf{Y}}\mathbf{Z}$$(20) where Z is a vector of independent random variables ${\mathrm{Z}}_i$ with $\mathbb{E}[{\mathrm{Z}}_i]=0$ and $\text{var}[{\mathrm{Z}}_i]=1$. Note, assigning a distributional form to the ${\mathrm{Z}}_i$ is a further choice; for example, we could assign a Gaussian, uniform, or any other distribution that would make sense given the context. Similar to how ensemble design is considered in history matching (e.g., Salter et al. 2019), we may eschew such distribution assumptions and set the bounds of the ${\mathrm{Z}}_i$ with an appeal to Chebyshev’s inequality. As is standard in the history matching literature, we set the concentration parameter k = 3 and sample ${\mathrm{Z}}_i\sim \mathcal{U}(-k,k)$; we call this the plausible set, and stress that it is not a probabilistic design, rather, a bounding notion of plausibility. To generate a joint sample (${\tilde{\mathbf{X}}}^{*},{\tilde{\mathbf{Y}}}^{*}$), we first sample ${\tilde{\mathbf{X}}}^{*}$ and then ${\tilde{\mathbf{Y}}}^{*}$. The structural dependencies between the sampled ${\tilde{\mathbf{X}}}^{*}$ and ${\tilde{\mathbf{Y}}}^{*}$ are captured by our updated beliefs on ${\mathcal{M}}_{\mathit{\beta }}$ and ${\mathbf{U}}_{\Theta }$. An example of a joint sample $({\tilde{\mathbf{X}}}^{*},{\tilde{\mathbf{Y}}}^{*})$ for the months of February and August is given in Figure 6.

Fig. 6 A single plausible sample from the SST and SIC reconstructions. (a) and (b) show SSTs from February and August, respectively and (c) and (d) show SICs from February and August, respectively.

We generate an ensemble of 25 runs comprising a single reference run using the mean SST and SIC fields produced as part of the PMIP4 LGM experiments, and 24 randomly generated plausible samples of SST and SIC. Averaging over PMIP models to act as boundary conditions in ice-sheet modeling is commonplace (e.g., Kageyama et al. 2017). Crucially, it should be noted that the reference run and each of the PMIP runs $(\mathbb{X},\mathbb{Y})$ do not lie within the calculated plausibility bounds. To examine the impact of running with plausible boundary conditions, we compare reference and plausible ice sheet heights at four geographically distinct locations: Arctic Canada, Central Greenland, Hudson Bay, and the Pacific coast; this is shown in Figure 7. The simulator is run for 5000 years, with the ice sheet initialized with the LGM Glac-1D reconstruction (Tarasov et al. 2012). Beyond the SST and SIC fields, no other model parameters were changed and the model set-up was based on previous simulations of the Greenland ice sheet (Gregory, George, and Smith 2020).

Fig. 7 Ensemble time-series plots at four spatial locations: Arctic Canada (A), Central Greenland (B), Hudson Bay (C), and the Pacific Coast (D). The reference run is shown in red, the simulations forced with sampled boundary conditions in black.

Figure 7 shows that the boundary conditions have a strong influence on the simulated ice sheet. SST and SIC can affect ice sheets either through changing the evaporation over the ocean that transforms into snow falling on the ice sheet, or by warming/cooling the regions close to the oceans thus affecting the surface melt rate. We find that the primary differences in ice sheet size are due to changes in snow accumulation. Our ensemble mainly produces lower ice elevation than the reference run. Examination of the individual ensemble members revealed that this was due to the cooler Eastern Pacific and Western Atlantic boundary conditions reducing evaporation over these oceans causing a reduction in the accumulation of snow onto the ice sheet. The difference between the reference run and the samples is most pronounced at the Pacific coast. This matches our expectation that ocean-proximal sites are more sensitive to marine influence of the SST and SIC than more continental sites. This is also where we see the strongest reduction in SST in our samples compared to the PMIP model simulations. Indeed, the second update causes a strong Pacific cooling along the coast of North America (Figure 2(e)) relative to what we expect based on the PMIP models. This is a region where models tend to underestimate the upwelling of cold waters from the deep ocean.

Such biases are common in climate models, but are particularly problematic for coupled climate-ice sheet models where the strong feedbacks between climate and ice sheets amplify the effects of climate biases, which can lead to runaway ice sheets (amplified growth or decay) and unrealistic geometries. The ensemble shows a substantial spread of ice elevation (5%–10% of height) caused by the variance in reconstructed SST and SIC, highlighting the importance for considering this source of uncertainty for modeling past ice sheets. Our results show that correcting for biases and incorporating uncertainty in surface ocean conditions has a substantial effect on the simulated ice sheet, which in turn influences the internal dynamics of the ice sheet and its vulnerability to collapse or propensity to grow. The ice sheet geometry itself has direct impact on atmospheric circulation and an indirect influence on ocean circulation from runoff, thus directly impacting global heat distributions and surface climate conditions. Our methodology provides a way to reduce climate biases by prescribing ocean surface conditions compatible with observations, while at the same time exploring the effects of this source of uncertainty on other parts of the earth system.

7 Discussion

By exploiting natural conditional exchangeability judgments we develop theory for the coexchangeable process model, as an extension to Rougier, Goldstein, and House (2013), that combines multi-model ensembles and data to model correlated spatio-temporal processes. We provide results for efficient and scalable inference that may be used for large spatio-temporal problems where probabilistic Bayesian methods are often computationally infeasible (see, e.g., Sansom, Stephenson, and Bracegirdle 2021). Our methodology requires fewer assumptions and less onerous belief specifications than that required by a probabilistic Bayesian analysis. To achieve these advancements we develop a Bayes linear analogue to a hierarchical Bayesian model. By combining exchangeability judgments and using the reparameterization of Hodges (1998), we extend the Bayes linear exchangeable regressions methodology. We obtain hitherto missing desirable properties present in traditional Bayesian hierarchical models such as the ability to make individual group level inference.

Large scale computational models often have complex spatio-temporal boundary conditions. This is particularly true for Earth system modeling, when any simulation of part of the Earth system, requires other spatio-temporal fields as boundary conditions. Our application looked at palaeo-era ice-sheet modeling, where our model had a coupled atmosphere and ice-sheet, with the SST and SIC as prescribed boundary conditions. These are usually specified using results from a reference simulation, or using a member of a Model Intercomparison Project (MIP). However, individual simulations of Earth system components are known to have biases and any individual simulation cannot adequately represent uncertainty due to boundary conditions. An idea for future investigation is to use the differences between MIP phases to estimate the ensemble discrepancy. This would allow for differences between older and newer MIP phases to inform the discrepancy between the newer models and reality. Whilst outside the scope of this work, using MIP iterations to inform model discrepancy is an interesting avenue, especially for models of the present-day where lots of data are available for validation.

Our methodology allows MIP simulations to be combined with observations efficiently, leading to joint reconstructions of climate boundary conditions that can be used in any area of Earth System modeling. We demonstrate its efficacy by reconstructing last glacial maximum SST and SIC to force an ice-sheet model. We show that the differences between reference ice-sheets and ice-sheets under plausible boundary conditions were considerable and that the uncertainty in the ice-sheet due to propagated boundary condition uncertainty is not ignorable. Other aspects of the Earth system are likely to be sensitive to their boundary conditions, so that joint reconstructions of the type we present here would allow MIP simulations and data to be combined in order to correct existing biases and quantify an important source of uncertainty. The use of MIP ensembles to drive simulations of Earth system components leads to important questions around how these ensembles should be designed. Our method makes the case that a priori exchangeability across as many models as possible is an important design goal.

Supplementary Materials

The supplementary material contains coexchangeable model and coexchangeable process model updating equations, calculations for belief adjustments, a proof to Theorem 4.1, and prior specifications for the application in Section 5.

Acknowledgments

Climate-ice sheet simulations were undertaken on ARC4, part of the High Performance Computing Facilities at the University of Leeds.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

All authors were funded by UKRI Future Leaders Fellowship MR/S016961/1.