# Coexchangeable process modelling for uncertainty quantification in joint
climate reconstruction

Lachlan Astfalck  School of Earth and Environment, The University of Leeds,
Leeds, UK  Daniel Williamson  Department of Mathematical Sciences, Exeter
University, Exeter, UK  The Alan Turing Institute, British Library, London, UK
Niall Gandy  School of Earth and Environment, The University of Leeds, Leeds,
UK  Lauren Gregoire  School of Earth and Environment, The University of Leeds,
Leeds, UK  Ruza Ivanovic  School of Earth and Environment, The University of
Leeds, Leeds, UK

###### 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. Whilst these boundary conditions are typically fixed using
available reconstructions in climate modelling studies, 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
coexchangable process model to capture multiple correlated spatio-temporal
fields of variables. We demonstrate that further exchangeability judgements
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 (19-23 ka) 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.

Keywords: Bayes linear methods, exchangeability analysis, multi model ensemble

##  1 Introduction

Numerical experiments are vital tools to climate science. Knowledge of
physical processes is embedded into software to simulate large scale events
that are otherwise impossible to observe at the required spatial and temporal
resolutions. The overwhelming majority of numerical climate simulators utilise
boundary conditions to represent non-calculated 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 historical
climate variables such as temperature, and precipitation (Gregoire et al.,
2012) that can be obtained using the Paleoclimate Model Intercomparison
Project (PMIP) ensemble of paleo-climate model runs (Ivanovic et al., 2016).
MMEs can be biased and do not necessarily span the uncertainty of the physical
process.

Some effort has been made to statistically reconstruct boundary conditions for
computer models. Salter et al. (2018) use 16 realisations from a single model
to develop a basis representation of a single spatio-temporal boundary
condition that could be sampled from and calibrated. However, in terms of the
fields that could be represented, this method was limited by the span of the
ensemble. Liu and Guillas (2017) capture a 3200 dimensional bathymetry
boundary condition for a Tsunami model using a dimension reduced kernel
approximation to an SPDE model that is trained using a large data set of ship
measurements. Whilst such approaches are attractive, they require a great deal
of data to accurately capture the boundary conditions, and these data are not
always available.

If modelling 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 et al. (2013) present a generalised framework,
therein termed the coexchangeable model, where exchangeability judgements 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 modelling task. The
coexchangeable model has seen some interest in more recent research, such as
in Sansom et al. (2021) where it is used to model emergent constraints for
future climate projections. Alternative to the methodology of Rougier et al.
(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 et al. (2013),
the model in Chandler (2013) is a special case of the coexchangeable model.

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 (we introduce this further in Section 4), where
the dependence between the boundary conditions is critically important. The
need for scalable inference and generality makes a coexchangeable approach
desirable, and thus we develop a hierarchical coexchangeable model that offers
a Bayes linear analogy to the natural Bayesian hierarchical model for the
problem. From simple and natural exchangeability judgements, we develop
efficient, scalable inference for joint reconstruction. This article proceeds
as follows. Section 2 reviews the existing statistical literature on Bayes
linear statistics and exchangeability analysis. Section 3 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 4 introduces our application to
reconstructing SST and SIC using our methodology. Section 5 discusses the
results of using the reconstructions of SST and SIC to simulate ice sheet-
atmosphere interactions at the LGM, and Section 6 concludes. Code and data are
provided as part of an R package downloadable from
github.com/astfalckl/exanalysis.

##  2 Review of statistical methodology

###  2.1 Bayes Linear Theory

Recent discourse increasingly acknowledges the difficulty of specifying
probability distributions to define likelihoods and prior beliefs for Bayesian
analyses (Barthelmé and Chopin, 2014; Bissiri et al., 2016; Jesus and
Chandler, 2011). This is especially relevant in climate applications where our
quantities of interest are large spatio-temporal fields for which data may be
(relatively) sparse. An alternative view treats expectation (rather than
probability) as the primitive quantity (De Finetti, 1975). Under expectation
as primitive, probabilities are the expectations of indicator functions for
events. This motivates second-order approaches such as Bayes linear methods
(see below and Goldstein and Wooff (2007)), where inference concerns
expectations and variances directly rather than as a by-product of
probabilistic inference.

Under the Bayes linear paradigm, beliefs on random quantities are described
via expectations and variances which are then adjusted by data. Consider
random quantity, $X$, with observations $X_{i}$, defined on the Hilbert space
$\mathcal{X}=\sum_{i}h_{i}X_{i}$ for
$X_{i}\in\mathcal{X}\subseteq\mathbb{R}^{n}$ and $h_{i}\in\mathbb{R}$, endowed
with inner product $\langle X,Y\rangle=\mathbb{E}[X^{\intercal}Y]$. Denote by
$D$ the collection of observed data, $D=(X_{1},\dots,X_{m})$, and subspace
$\mathcal{D}$ as $\mathcal{D}=\text{span}[1,D]$. The adjusted expectation
(written in terms of $X$ and $D$ as $\mathbb{E}_{D}[X]$, that is, the
expectation of our beliefs $X$ adjusted by data $D$) is the element in
$\mathcal{D}$ that minimises $\|X-\mathbb{E}_{D}[X]\|$. The adjusted
expectation has solution

$\mathbb{E}_{D}[X]=\mathbb{E}[X]+\text{cov}[X,D]\text{var}[D]^{\dagger}(D-\mathbb{E}[D])$
(1)

where $\text{var}[D]^{\dagger}$ is any pseudo-inverse of $\text{var}[D]$, most
commonly the Moore-Penrose inverse. Equation (1) describes the orthogonal
projection of each element of our beliefs, $X$, onto $\mathcal{D}$. The
adjusted variance, $\text{var}_{D}[X]$, is the outer product
$\mathbb{E}\left[\left(X-\mathbb{E}_{D}[X]\right)\left(X-\mathbb{E}_{D}[X]\right)^{\intercal}\right]$
given $\mathbb{E}_{D}[X]$ in (1) and is

$\text{var}_{D}[X]=\text{var}[X]-\text{cov}[X,D]\text{var}[D]^{\dagger}\text{cov}[D,X].$
(2)

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

####  2.1.1 Bayes linear analysis of exchangeable data

Exchangeability for a sequence of random quantities within a probabilistic
Bayesian analysis represents a simple symmetry judgement that implies that any
finite sub-collection of quantities within the sequence have the same
distribution (a priori). Second-order exchangeability for such a sequence is a
similar analogue for judgements when expectation is primitive, and implies
that each element of the same sequence has the same (prior) expectation and
variance and that all pairs have the same covariance. If data, $D$, are
second-order exchangeable we may, according to the second-order representation
theorem (Goldstein, 1986), write

$X_{i}=\mathcal{M}(X)+\mathcal{R}_{i}(X)$ (3)

for a common mean $\mathcal{M}(X)$ and uncorrelated residuals
$\mathcal{R}_{i}(X)$. Second-order exchangeability implies
$\text{cov}[X_{i},X_{i^{\prime}}]=\Psi+\mathbf{1}_{i=i^{\prime}}\Gamma\;\forall
i,i^{\prime}$. For second-order exchangeable sequences there is predictive
sufficiency for updating beliefs on $X$ by only updating $\mathcal{M}(X)$ by
the data $D$. Mathematically we write this as
$[D\;\rotatebox[origin={c}]{90.0}{$\models$}\;X]/\mathcal{M}(X)$, i.e., given
$\mathcal{M}(X)$, $D$ and $X$ are uncorrelated. Further, the sample mean,
$\bar{X}=\frac{1}{m}\sum_{i=1}^{m}X_{i}$, is Bayes linear sufficient for
updating beliefs on $\mathcal{M}(X)$ (and consequently on $X$):
$[D\;\rotatebox[origin={c}]{90.0}{$\models$}\;\mathcal{M}(X)]/\bar{X}$. A
proof of Bayes linear sufficiency for exchangeable data is provided in Section
6.10 of Goldstein and Wooff (2007).

###  2.2 Exchangeability analysis for multi-model ensembles

We now review the coexchangeable model of Rougier et al. (2013) who use a
Bayes linear approach to describe MMEs as follows. Let
$\mathbb{X}\coloneqq\\{\textbf{X}_{1},\dots,\textbf{X}_{m}\\}$ be a collection
of $q$-dimensional outputs from the $m$ simulators that form the MME;
$\textbf{X}^{*}$, the true unobserved process that the simulators aim to
model; and $\textbf{Z}_{\textbf{X}}$, the noisy and incomplete measurement of
$\textbf{X}^{*}$ in the real world. The model requires only two assumptions:
first, that the $\textbf{X}_{i}$ are second-order exchangeable and, second,
that the $\textbf{X}_{i}$ are coexchangeable with the truth $\textbf{X}^{*}$,
implying $\text{cov}[\textbf{X}^{*},\textbf{X}_{i}]=\Sigma,\;\forall i$.

As in (3), exchangeability within $\mathbb{X}$ implies

$\textbf{X}_{i}=\mathcal{M}(\textbf{X})+\mathcal{R}_{i}(\textbf{X}),\hskip
14.22636pti=1,\dots,m,$ (4)

where $\mathcal{M}(\textbf{X})$ is a shared mean term, the
$\mathcal{R}_{i}(\textbf{X})$ are the zero-mean, uncorrelated residuals of
each simulator with common variance, and the $\mathcal{M}(\textbf{X})$ and
$\mathcal{R}_{i}(\textbf{X})$ are uncorrelated. Coexchangeability between
$\textbf{X}^{*}$ and $\mathbb{X}$ implies sufficiency of
$\mathcal{M}(\textbf{X})$ for $\textbf{X}^{*}$,
$[\textbf{X}_{i}\;\rotatebox[origin={c}]{90.0}{$\models$}\;\textbf{X}^{*}]/\mathcal{M}(\textbf{X})\;\forall
i$. This allows us to write

$\textbf{X}^{*}=\textbf{A}\mathcal{M}(\textbf{X})+\textbf{U}_{\textbf{X}}$ (5)

where A is a known matrix, and $\textbf{U}_{\textbf{X}}$ represents the
ensemble discrepancy that is uncorrelated with $\mathcal{M}(\textbf{X})$ and
the $\mathcal{R}_{i}(\textbf{X})$. The data, $\textbf{Z}_{\textbf{X}}$ are
modelled as

$\textbf{Z}_{\textbf{X}}=\textbf{H}_{\textbf{X}}\textbf{X}^{*}+\textbf{W}_{\textbf{X}}$
(6)

for measurement error $\textbf{W}_{\textbf{X}}$, and known incidence matrix
$\textbf{H}_{\textbf{X}}$. The statistical model defined by (4)–(6) will
hereafter by referred to as the “coexchangeable model”. Inference for this
model makes use of Bayes linear sufficiency of the ensemble mean,
$\bar{\textbf{X}}=\frac{1}{m}\sum_{i=1}^{m}\textbf{X}_{i}$, for updating by
$\mathbb{X}$ and is therefore very efficient. We outline the inferential
procedure developed by Rougier et al. (2013) in the supplementary material.

##  3 The coexchangeable process model

Now consider each simulator as producing output pairs of $q$-dimensional
fields $(\textbf{X}_{i},\textbf{Y}_{i})$ for $i=1,\ldots,m$, with the
corresponding fields in reality denoted $(\textbf{X}^{*},\textbf{Y}^{*})$, for
which we have (partial) observations $\textbf{Z}_{\textbf{X}}$ and
$\textbf{Z}_{\textbf{Y}}$, made with some error. Whilst it is possible to
consider the coexchangeable model over both fields simultaneously, a
probabilistic Bayesian approach would normally leverage conditional
probability laws to capture structured dependencies between the fields (these
are a direct consequence of certain exchangeability judgements as with Sansom
et al. (2021)). In what follows we capture explicit structural dependencies
between random quantities $(\textbf{X},\textbf{Y})$, by first only imposing
the coexchangeable model for $\textbf{X}_{i},\textbf{X}^{*}$, and then
considering exchangeability judgements over the processes
$\textbf{Y}_{i}(\textbf{X}_{i})$ and coexchangeability of
$\textbf{Y}^{*}(\textbf{X}^{*})$.

The rest of this section develops theory and methodology to model the random
process $\textbf{Y}(\textbf{X})$ within a coexchangeability framework. Section
3.1 develops a general theory for second order exchangeable processes as a
Bayes linear analogy to a Bayesian hierarchical model. Section 3.2 exploits
the underlying geometry induced by the exchangeability for efficient updating.
Finally, Section 3.3 develops our coexchangeable process model for the stated
problem and gives results for efficient inference.

###  3.1 Exchangeable processes

Assume we observe a sequence of observation/covariate pairings from group $i$,
$(Y_{i,1},X_{i,1}),(Y_{i,2},X_{i,2}),\dots$. Here we consider each simulator
in our MME to comprise one such group and for the $(Y_{i,j},X_{i,j})$ to be
some $p$-dimensional sub-level process contained in
$(\textbf{Y}_{i},\textbf{X}_{i})$. To help fix ideas, in the application in
Section 4, $\textbf{Y}_{i}$ is spatio-temporal SIC and $\textbf{X}_{i}$ is
spatio-temporal SST, both for simulator $i$. Each $j$ indexes time (so that
the $Y_{i,j}$ and $X_{i,j}$ are spatial processes), though it could instead
index space, both space and time, or some other feature of the process. We
impose conditional second order exchangeability within each group, the second
order equivalent of the exchangeability judgements used in Bayesian regression
(see Williamson and Sansom, 2019, for discussion). Specifically, conditional
second order exchangeability implies that the $Y_{i,1},Y_{i,2},...$ are second
order exchangeable given fixed $X_{i,j}=X_{i}\;\forall j$. We may thus specify
a mean process,
$\mathbb{E}[Y_{i,j}(X_{i,j})]=\mathcal{M}(X_{i,j},\beta_{i})\in\mathbb{R}^{p}$
for some known covariate $X_{i,j}$ and unknown parameter
$\beta_{i}\in\mathbb{R}^{k}$ specific to group $i$. According to the second
order representation theorem,

$Y_{i,j}(X_{i,j})=\mathcal{M}(X_{i,j},\beta_{i})+\mathcal{R}_{i,j}(Y_{i,j};X_{i,j})$
(7)

and
$\text{cov}[Y_{i,j}(X_{i,j}),Y_{i,j^{\prime}}(X_{i,j^{\prime}})]=\text{cov}[\mathcal{M}(X_{i,j},\beta_{i}),\mathcal{M}(X_{i,j^{\prime}},\beta_{i})]+\mathbf{1}_{\\{j=j^{\prime}\\}}\text{var}[\mathcal{R}_{i,j}(Y_{i,j};X_{i,j})]$.
$\mathcal{M}(X_{i,j},\beta_{i})$ may represent any functional relationship
between $X_{i,j}$ and $\beta_{i}$ and, to infer $\beta_{i}$ from
$Y_{i,j}(X_{i,j})$, we need only specify a joint inner product space in which
they reside. We set

$Y_{i,j}(X_{i,j})=\phi(X_{i,j})\beta_{i}+\mathcal{R}_{i,j}(Y_{i,j};X_{i,j})$
(8)

where $\phi(X_{i,j})$ is a matrix with $k$ columns representing known basis
functions of $X_{i,j}$. The choice to model $\mathcal{M}(X_{i,j},\beta_{i})$
linear in $\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 the
random vectors $\beta_{1},\beta_{2},\dots$ so that

$\beta_{i}=\mathcal{M}(\beta)+\mathcal{R}_{i}(\beta),$ (9)

with expectation $\mathbb{E}[\beta_{i}]=\mathcal{M}(\beta)$ and covariance
$\text{cov}[\beta_{i},\beta_{i^{\prime}}]=\text{cov}[\mathcal{M}(\beta_{i}),\mathcal{M}(\beta_{i^{\prime}})]+\mathbf{1}_{\\{i=i^{\prime}\\}}\text{var}[\mathcal{R}_{i}(\beta)]$,
by the representation theorem.

Observe that
$\textbf{Y}_{i}=(Y_{i,1}(X_{i,1})^{\intercal},\dots,Y_{i,n}(X_{i,n})^{\intercal})^{\intercal}$
and
$\textbf{X}_{i}=(X_{i,1}^{\intercal},\dots,X_{i,n}^{\intercal})^{\intercal}$
for $j\in\\{1,\dots,n\\}$, and define $\Phi_{i}\coloneqq\phi(\textbf{X}_{i})$,
and
$\mathcal{R}_{i}(\textbf{Y}_{i};\textbf{X}_{i})\coloneqq(\mathcal{R}_{i,1}(Y_{i,1};X_{i,1})^{\intercal},\dots,\mathcal{R}_{i,n}(Y_{i,n};X_{i,n})^{\intercal})^{\intercal}$.
We make additional natural belief specifications,
$\text{cov}[\mathcal{R}_{i}(\textbf{Y}_{i};\textbf{X}_{i}),\mathcal{R}_{i^{\prime}}(\beta)]=0$,
and
$\text{cov}[\mathcal{R}_{i}(\textbf{Y}_{i};\textbf{X}_{i}),\mathcal{M}(\beta)]=0\;\forall
i,i^{\prime}$, and further note that
$\text{cov}[Y_{i,j}(X_{i,j}),\beta_{i^{\prime}}]=\phi(X_{i,j})\text{cov}[\beta_{i},\beta_{i^{\prime}}]$
results from our exchangeability judgements. Equations (8) and (9) imply

$\small\normalsize\begin{bmatrix}\textbf{Y}_{1}\\\ \vdots\\\ \textbf{Y}_{m}\\\
\beta_{1}\\\ \vdots\\\
\beta_{m}\end{bmatrix}=\left[\begin{array}[]{@{}c|c@{}}\begin{matrix}\Phi_{1}&\cdots&\bm{0}_{q\times
k}\\\ \vdots&\ddots&\vdots\\\ \bm{0}_{q\times
k}&\cdots&\Phi_{m}\end{matrix}&\bm{0}_{qm\times k}\\\
\hline\cr\bm{0}_{km\times km}&\begin{matrix}\textbf{J}_{m\times
1}\otimes\textbf{I}_{k}\end{matrix}\end{array}\right]\begin{bmatrix}\beta_{1}\\\
\vdots\\\ \beta_{m}\\\
\mathcal{M}(\beta)\end{bmatrix}+\begin{bmatrix}\mathcal{R}_{1}(\textbf{Y}_{1};\textbf{X}_{1})\\\
\vdots\\\ \mathcal{R}_{m}(\textbf{Y}_{m};\textbf{X}_{m})\\\
\mathcal{R}_{1}(\beta)\\\ \vdots\\\ \mathcal{R}_{m}(\beta)\end{bmatrix},$ (10)

where $\bm{0}_{a\times b}$ is a $a\times b$ zero matrix, $\textbf{J}_{a\times
b}$ is a $a\times b$ matrix of ones, and $\textbf{I}_{a}$ is the $a\times a$
identity matrix.

We note that the specification of (10) differs from the standard approach to
modelling exchangeable processes in Bayes linear statistics where inference
only concerns $\mathcal{M}(\beta)$ by, in effect, substituting (9) into (8)
(see Goldstein and Wooff (1998)). By jointly modelling the $\beta_{i}$ and
$\mathcal{M}(\beta)$ we provide a closer analogy to Bayesian hierarchical
models. Advantages of (10) are that group level inference on the $\beta_{i}$
are available, and our adjusted beliefs are subject to the same types of
shrinkage effects seen in Bayesian hierarchical models (for instance, see
Section 8.4.2 of Hoff (2009)). In the context of climate modelling, revised
beliefs of the $\beta_{i}$ provide inference on the physical processes at the
model level.

###  3.2 Inference for exchangeable processes

Noting that (9) is equivalently restated as
$0=\mathcal{M}(\beta)+\mathcal{R}_{j}(\beta)-\beta_{j}$, and using the re-
parametrisation of Hodges (1998), we can re-express (10) via

$\small\normalsize\begin{bmatrix}\textbf{Y}_{1}\\\ \vdots\\\ \textbf{Y}_{m}\\\
\bm{0}_{km\times
1}\end{bmatrix}=\left[\begin{array}[]{@{}c|c@{}}\begin{matrix}\Phi_{1}&\cdots&\bm{0}_{q\times
k}\\\ \vdots&\ddots&\vdots\\\ \bm{0}_{q\times
k}&\cdots&\Phi_{m}\end{matrix}&\bm{0}_{qm\times k}\\\
\hline\cr-\textbf{I}_{km}&\begin{matrix}\textbf{J}_{m\times
1}\otimes\textbf{I}_{k}\end{matrix}\end{array}\right]\begin{bmatrix}\beta_{1}\\\
\vdots\\\ \beta_{m}\\\
\mathcal{M}(\beta)\end{bmatrix}+\begin{bmatrix}\mathcal{R}_{1}(\textbf{Y}_{1};\textbf{X}_{1})\\\
\vdots\\\ \mathcal{R}_{m}(\textbf{Y}_{m};\textbf{X}_{m})\\\
\mathcal{R}_{1}(\beta)\\\ \vdots\\\ \mathcal{R}_{m}(\beta)\end{bmatrix}$ (11)

with the familiar linear canonical form
$\textbf{Y}=\textbf{X}\textbf{B}+\textbf{E}$. As we consider the joint
representation, expectation is taken jointly over $\textbf{Y}_{i}$ and
$\beta_{i}$. Belief specifications follow from Section 3.1 and the adjusted
expectation, $\mathbb{E}_{\textbf{Y}}[\textbf{B}]$, and adjusted variance,
$\text{var}_{\textbf{Y}}[\textbf{B}]$, are calculated via (1) and (2); the
joint belief specifications and specific updating equations are provided in
the supplementary material.

Calculating $\mathbb{E}_{\textbf{Y}}[\textbf{B}]$ and
$\text{var}_{\textbf{Y}}[\textbf{B}]$ requires the inversion of the
$m(q+k)\times m(q+k)$ matrix, $\text{var}[\textbf{Y}]$. We show in Theorem 3.1
(below) that a Bayes linear sufficiency argument can be made to make
equivalent inference for B instead with the inversion of a smaller (when
$k<q$) $2km\times 2km$ matrix. For climate models, $q$ is generally large as
either the size of the process, $p$, is large (such as when $Y_{i,j}$ is a
high-resolution spatial process), or $n$ is large (such as when $j$ indexes
space and/or time). For many basis designs $k<<q$, and when groups, $i$, index
MIP simulations, $m$ is generally small (our application has $m=13$).

######  Theorem 3.1.

Let $\hat{\bm{\beta}}=(\hat{\beta}_{1},\dots,\hat{\beta}_{m})$ with
$\hat{\beta}_{i}=(\Phi_{i}^{\intercal}\Phi_{i})^{\dagger}\Phi_{i}^{T}\textbf{Y}_{i}$.
Then $\hat{\bm{\beta}}$ is Bayes linear sufficient for Y for adjusting B if
the column space of projection matrix
$P_{i}=\Phi_{i}(\Phi_{i}^{\intercal}\Phi_{i})^{\dagger}\Phi_{i}^{\intercal}$,
$C(P_{i})$, is invariant over $i$, i.e., $C(P_{i})=\mathcal{C}$ $\forall i$.

###### Proof.

As $P_{i}$ is a projection matrix (this is immediately verified as
$P_{i}^{2}=P_{i}$), $\hat{\textbf{Y}}_{i}=P_{i}\textbf{Y}_{i}\in C(P_{i})$
projects $\textbf{Y}_{i}$ into the column space of $P_{i}$ and
$\hat{R}_{i}=\textbf{Y}_{i}-\hat{\textbf{Y}}_{i}=(1-P_{i})\textbf{Y}_{i}\in
N(P_{i})$ projects $\textbf{Y}_{i}$ into the null space of $P_{i}$,
$N(P_{i})$. Define $\mathcal{Y}_{i}$ as the inner product space in which
$\textbf{Y}_{i}$ resides and on which $P_{i}$ is defined. $\mathcal{Y}_{i}$
may be partitioned as direct sum $\mathcal{Y}_{i}=C(P_{i})\oplus N(P_{i})$. As
$C(P_{i})$ and $N(P_{i})$ form orthogonal subspaces, $\langle
A,B\rangle_{\mathcal{Y}_{i}}=0$ and so $\text{cov}[A,B]=0$ for any vector
$A\in C(P_{i})$ and $B\in N(P_{i})$.

From (11), the conditions
$[\textbf{Y}_{i}\;\rotatebox[origin={c}]{90.0}{$\models$}\;\mathcal{M}(\beta)]/\hat{\beta}_{i}$
and
$[\textbf{Y}_{i}\;\rotatebox[origin={c}]{90.0}{$\models$}\;\beta_{i^{\prime}}]/\hat{\beta}_{i}\;\forall
i,i^{\prime}$ satisfy Bayes linear sufficiency,
$[\textbf{Y}\;\rotatebox[origin={c}]{90.0}{$\models$}\;\textbf{B}]/\hat{\bm{\beta}}$.
From Theorem 5.20.7 Goldstein and Wooff (2007),
$[\textbf{Y}_{i}\;\rotatebox[origin={c}]{90.0}{$\models$}\;\mathcal{M}(\beta)]/\hat{\beta}_{i}$
and
$[\textbf{Y}_{i}\;\rotatebox[origin={c}]{90.0}{$\models$}\;\beta_{i^{\prime}}]/\hat{\beta}_{i}\;\forall
i,i^{\prime}$ are satisfied if
$\text{cov}[\textbf{Y}_{i}-\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}],\mathcal{M}(\beta)]=0$
and
$\text{cov}[\textbf{Y}_{i}-\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}],\beta_{i^{\prime}}]=0$,
respectively. Observe that
$\text{cov}[\textbf{Y}_{i},\hat{\textbf{Y}}_{i^{\prime}}]=\text{cov}[P_{i}\textbf{Y}_{i}+(1-P_{i})\textbf{Y}_{i},\hat{\textbf{Y}}_{i^{\prime}}]=\text{cov}[\hat{\textbf{Y}}_{i},\hat{\textbf{Y}}_{i^{\prime}}]\;\forall
i,i^{\prime}$ so that
$\text{cov}[\textbf{Y}_{i},\hat{\bm{\beta}}_{i^{\prime}}]=\Phi_{i}\text{cov}[\hat{\bm{\beta}}_{i},\hat{\bm{\beta}}_{i^{\prime}}]\;\forall
i,i^{\prime}$ and thus calculating
$\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}]$ via (1) yields
$\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}]=\Phi_{i}\hat{\beta}_{i}=\hat{\textbf{Y}}_{i}$.
If $C(P_{i})$ and $N(P_{i})$ are invariant over $i$ (i.e.,
$C(P_{i})=\mathcal{C}$ and $N(P_{i})=\mathcal{N}$) then $C(P_{i})\perp
N(P_{i^{\prime}})$ $\forall i,i^{\prime}$. Thus
$\textbf{Y}_{i}-\hat{\textbf{Y}}_{i}\in\mathcal{N}$ and
$\Phi_{i}\beta_{i},\Phi_{i}\mathcal{M}(\beta)\in\mathcal{C},\;\forall i$, so
that
$\text{cov}[\textbf{Y}_{i}-\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}],\Phi_{i}\mathcal{M}(\beta)]=0$
and
$\text{cov}[\textbf{Y}_{i}-\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}],\Phi_{i}\beta_{i^{\prime}}]=0\;\forall
i,i^{\prime}$, and thus
$\text{cov}[\textbf{Y}_{i}-\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}],\mathcal{M}(\beta)]=0$
and
$\text{cov}[\textbf{Y}_{i}-\mathbb{E}_{\hat{\beta}_{i}}[\textbf{Y}_{i}],\beta_{i^{\prime}}]=0\;\forall
i,i^{\prime}$. ∎

As $C(P_{i})=C(\Phi_{i})$, if $\Phi_{i}$ is full rank this is sufficient (but
not necessary) in satisfying the condition $C(P_{i})=\mathcal{C}\;\forall i$
in Theorem 3.1 — most sensible basis designs will ensure this. Assuming
$C(P_{i})=\mathcal{C}\;\forall i$, and appealing to the sufficiency of
$\hat{\bm{\beta}}$ for Y for adjusting B, we write the joint specification of
the exchangeability judgements made in Section 3.1 as

$\small\normalsize\begin{bmatrix}\hat{\beta}_{1}\\\ \vdots\\\
\hat{\beta}_{m}\\\ \bm{0}_{km\times
1}\end{bmatrix}=\left[\begin{array}[]{@{}c|c@{}}\begin{matrix}\textbf{I}_{km}\end{matrix}&\bm{0}_{km\times
k}\\\ \hline\cr-\textbf{I}_{km}&\begin{matrix}\textbf{J}_{m\times
1}\otimes\textbf{I}_{k}\end{matrix}\end{array}\right]\begin{bmatrix}\beta_{1}\\\
\vdots\\\ \beta_{m}\\\
\mathcal{M}(\beta)\end{bmatrix}+\begin{bmatrix}\mathcal{R}_{1}(\hat{\beta};\textbf{X}_{1})\\\
\vdots\\\ \mathcal{R}_{m}(\hat{\beta};\textbf{X}_{m})\\\
\mathcal{R}_{1}(\beta)\\\ \vdots\\\ \mathcal{R}_{m}(\beta)\end{bmatrix}$ (12)

where
$\mathcal{R}_{i}(\hat{\beta};\textbf{X}_{i})=(\Phi_{i}^{\intercal}\Phi_{i})^{\dagger}\Phi_{i}^{\intercal}\mathcal{R}_{i}(\textbf{Y}_{i};\textbf{X}_{i})$.
Calculation of $\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{B}]$ and
$\text{var}_{\hat{\bm{\beta}}}[\textbf{B}]$ follows (1) and (2); the specific
equations are given in the supplementary material. If interest lies only in
the belief updates of $\mathcal{M}(\beta)$, then computation of
$\mathbb{E}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ and
$\text{var}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ exist by only computing
the requisite parts of $\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{B}]$ and
$\text{var}_{\hat{\bm{\beta}}}[\textbf{B}]$; equations for
$\mathbb{E}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ and
$\text{var}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ are also provided in the
supplementary material.

###  3.3 Reality as a coexchangeable process

We now show how judgements of coexchangeability may be made to incorporate
(11) or (12) into the coexchangeable model. Define the true process of
$\textbf{Y}_{i}$ that the simulators attempt to resolve given some
$\textbf{X}^{*}$ as $\textbf{Y}^{*}(\textbf{X}^{*})$. We impose
coexchangeability of $\textbf{Y}^{*}(\textbf{X}^{*})$ and $\textbf{Y}_{i}$
given fixed $\textbf{X}_{i}=\textbf{X}$ (and hence $\Phi_{i}=\Phi$) $\forall
i$, i.e.,
$\text{cov}[\textbf{Y}^{*}(\textbf{X}^{*}),\textbf{Y}_{i}]=\Sigma(\textbf{X}^{*})\;\forall
i$. As
$\text{cov}[\textbf{Y}^{*}(\textbf{X}^{*}),\textbf{Y}_{i}]=\text{cov}[\textbf{Y}^{*}(\textbf{X}^{*}),\Phi\beta_{i}]=\Sigma(\textbf{X}^{*})=\Gamma(\textbf{X}^{*})\Phi^{\intercal}$
we find that
$\text{cov}[\textbf{Y}^{*}(\textbf{X}^{*}),\beta_{i}]=\Gamma(\textbf{X}^{*})\;\forall
i$ and so we have coexchangeability of $\textbf{Y}^{*}(\textbf{X}^{*})$ and
$\beta_{i}$. Bayes linear sufficiency of $\mathcal{M}(\beta)$ for the
$\textbf{Y}_{i}$ for adjusting $\textbf{Y}^{*}(\textbf{X}^{*})$ follows. We
write a model for $\textbf{Y}^{*}(\textbf{X}^{*})$ such that

$\textbf{Y}^{*}(\textbf{X}^{*})=A(\textbf{X}^{*})\mathcal{M}(\beta)+\textbf{U}_{\textbf{Y}},$
(13)

where $\textbf{U}_{\textbf{Y}}$ is a model-mismatch term orthogonal to the
$\textbf{Y}_{i}$ and the assumption of coexchangeability permits
$A(\textbf{X}^{*})$ to be any matrix of suitable dimensions. The choice of
$A(\textbf{X}^{*})=\phi(\textbf{X}^{*})$ is logical in our context as an
analogue to (8) for the reality model, but different choices are permissible
should they make sense to the application. Finally, the data,
$\textbf{Z}_{\textbf{Y}}$ are modelled as

$\textbf{Z}_{\textbf{Y}}=\textbf{H}_{\textbf{Y}}\textbf{Y}^{*}(\textbf{X}^{*})+\textbf{W}_{\textbf{Y}},$
(14)

for measurement error $\textbf{W}_{\textbf{Y}}$, and known incidence matrix
$\textbf{H}_{\textbf{Y}}$.

Similar to Rougier et al. (2013), we assume sufficiency of
$\textbf{Y}^{*}(\textbf{X}^{*})$ for $\textbf{Z}_{\textbf{Y}}$, allowing us to
update beliefs on $\textbf{Y}^{*}(\textbf{X}^{*})$ in two stages. As adjusted
expectation is a linear operator (Property 3.2.1 Goldstein and Wooff (2007)),
$\textbf{U}_{\textbf{Y}}$ is defined uncorrelated with the ensemble, and
supposing the conditions of Theorem 3.1 hold, the adjustments of
$\textbf{Y}^{*}(\textbf{X}^{*})$ by the $(\textbf{Y}_{i},\textbf{X}_{i})$ are

$\displaystyle\mathbb{E}_{(\textbf{Y},\textbf{X})}[\textbf{Y}^{*}(\textbf{X}^{*})]$
$\displaystyle=\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]=\phi(\textbf{X}^{*})\mathbb{E}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]+\mathbb{E}[\textbf{U}_{\textbf{Y}}]\text{,
and}$ (15)
$\displaystyle\text{var}_{(\textbf{Y},\textbf{X})}[\textbf{Y}^{*}(\textbf{X}^{*})]$
$\displaystyle=\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]=\phi(\textbf{X}^{*})\text{var}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]\phi(\textbf{X}^{*})^{\intercal}+\text{var}[\textbf{U}_{\textbf{Y}}].$
(16)

Note that if the conditions of Theorem 3.1 are not met, then the first
equality is not necessarily true and the adjustments on the right hand side of
these equations would be with respect to $(\textbf{Y},\textbf{X})$. The second
update of $\textbf{Y}^{*}(\textbf{X}^{*})$ by data $\textbf{Z}_{\textbf{Y}}$
is calculated as

$\displaystyle\begin{split}\mathbb{E}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\textbf{X}^{*})]&=\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]+\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]\textbf{H}_{\textbf{Y}}^{\intercal}\\\
(\textbf{H}_{\textbf{Y}}\text{var}_{\hat{\bm{\beta}}}&[\textbf{Y}^{*}(\textbf{X}^{*})]\textbf{H}_{\textbf{Y}}^{\intercal}+\text{var}[\textbf{W}_{\textbf{Y}}])^{-1}(\textbf{Z}_{\textbf{Y}}-\textbf{H}_{\textbf{Y}}\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]),\text{
and}\end{split}$ (17)
$\displaystyle\begin{split}\text{var}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\textbf{X}^{*})]&=\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]-\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]\textbf{H}_{\textbf{Y}}^{\intercal}\\\
(\textbf{H}_{\textbf{Y}}\text{var}_{\hat{\bm{\beta}}}&[\textbf{Y}^{*}(\textbf{X}^{*})]\textbf{H}_{\textbf{Y}}^{\intercal}+\text{var}[\textbf{W}_{\textbf{Y}}])^{-1}\textbf{H}_{\textbf{Y}}\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})].\end{split}$
(18)

##  4 Joint reconstruction of Paleo sea-surface temperature and sea-ice with
PMIP3 and PMIP4

Our application is the joint reconstruction of SST and SIC in order to provide
spatio-temporal boundary conditions for coupled ice sheet-atmosphere
simulations. These simulations will be run with the FAMOUS-Ice simulator
(Smith et al., 2021), a computationally efficient model which can represent
ice sheet-climate processes in sufficient detail to accurately simulate the
contemporary Greenland ice sheet (Gregory et al., 2020). One way to test the
models’ ability to simulate climate and ice sheets different from today (e.g.
the future), is to simulate the past (Kageyama et al., 2017; Harrison et al.,
2015; Schmidt et al., 2014). The deglaciation ($\sim$21,000 to 7,000 years
ago) 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), which respectively warmed and shrunk due to increasing greenhouse gas
concentrations. To begin to study the deglaciation, an ice sheet needs to be
“grown” within the model under the steady state boundary conditions (SST and
SIC) of the relatively stable conditions of the LGM ($\sim$23,000 to 19,000
years ago). The character of the resulting ice sheets are sensitive to these
boundary conditions and so it is cruicial 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 include the Climate Long range
Investigation, Mapping, and Prediction (CLIMAP) project (CLIMAP Project,
1981), the Glacial Atlantic Ocean Mapping (GLAMAP) project (Sarnthein et al.,
2003), the Multiproxy Approach for the Reconstruction of the Glacial Ocean
Surface (MARGO) project (Kucera et al., 2005), and more recently Paul et al.
(2020). These global reconstructions (as well as other regionalised
reconstructions such as De Vernal et al. (2005) and Gersonde et al. (2005))
are ‘data only’ based solely on paleodata syntheses and do not include
numerical simulations in their analyses. As such, they are subject to large
measurement error, large biases in polar regions, incomplete spatial coverage,
and poor temporal resolution. Modern statistical methods (e.g. in Zhang et al.
(2020) and Sha et al. (2019)) for reconstructing spatio-temporal climate
fields rely on computationally intensive methods that require large
observational data and so are inappropriate here given the relative data
sparsity at the LGM.

Coupled simulations of the LGM are useful to ensure spatio-temporal
consistency of SST and SIC, but these can also be very different from
observations in critical regions for growing ice sheets (Salter et al., 2019).
Tierney et al. (2020) use a variant of the ensemble square-root Kalman filter
(Whitaker and Hamill, 2002) to adjust a model with observations to provide SST
reconstruction (with uncertainty quantification), based on the iCESM model.
Differences in model physics are typically much wider across models than the
variation of a particular model when varying its parameters (Masson and
Knutti, 2011). Therefore, reconstructions based on a single model may be
biased and overconfident. PMIP defines experimental protocols for coupled
ocean-atmosphere models (with prescribed greenhouse gases and ice sheets) to
simulate the LGM climate (Kageyama et al., 2017, 2021). The PMIP community has
produced MMEs for different phases of the project; PMIP2, PMIP3, PMIP4 run
with different generations of model and slight adjustments in input greenhouse
gases and ice sheets. Annan and Hargreaves (2013) made use of such MMEs to
produce an SST reconstruction by regressing the proxy data to the PMIP2 model
data. Other uses of PMIP data have been to directly force ice sheet
simulations (Gregoire et al., 2016).

Our methodology allows for the first joint reconstruction of SST and SIC (we
are aware of) that coherently combines the PMIP models (we use PMIP3 and
PMIP4) with available proxy data to deliver boundary conditions with
uncertainty quantification. We also demonstrate inference for the differences
in model physics that our Bayes linear analogy to the hierarchical model
offers.

###  4.1 Numerical simulators and paleo-observations

####  4.1.1 The MME: PMIP global climate models

As with other exchangeable MME methods (Rougier et al., 2013; Sansom et al.,
2021), we select a single representative simulation from each modelling group
that contributed to the PMIP3 and PMIP4 MMEs (Kageyama et al., 2021), in order
to make the assumption of prior exchangeability, and exchangeable processes,
reasonable. However, for the HadCM3 model simulations, we make use of all
three available PMIP4 simulations that use different ice sheet boundary
conditions as these differences are known to affect the simulated SST and SIC.
The $m=13$ models selected were

$\begin{split}S=\\{&\texttt{CNRM-3, IPSL-3, MIROC-3, MPI-3, CCSM4-3, GISS-3,
MRI-3, AWI-4,}\\\ &\texttt{MIROC-4, MPI-4, HadCM3-PMIP, HadCM3-Glac1D,
HadCM3-Ice6G}\\}.\end{split}$ (19)

Each ensemble member was projected onto the FAMOUS ocean grid with $p=4145$
spatial locations, and we use the 12 monthly means of SST and SIC. The choice
to use monthly means rather than a full simulation is based on only having
monthly means for $7$ of the $13$ ensemble members.

####  4.1.2 SST Observations

We use proxy temperature data, $\textbf{Z}_{\textbf{X}}$, obtained from sea-
bed sediment core samples to infer LGM SST, $\textbf{X}^{*}$. We use the MARGO
compilation (Kucera et al., 2005) supplemented with some more recent data from
Benz et al. (2016) to fill in some of the sparsity in the Southern Ocean. Most
proxy temperature data measure SST as an annual mean, but Southern Ocean
locations are provided as summer means. The measurements from final data
compilation are shown in Figure 1, annual and summer means are depicted with
points and triangles, respectively.

Figure 1: Proxy-based measurements of LGM SST and maximum sea-ice extents.
Annual mean and summer mean sea-surface temperatures 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.

There is very likely some strong observation bias in the foraminifera-based
proxy measurements from the Arctic and Nordic seas, approximately north of
$62^{\circ}$ N. For example, in this region, very cold water and full sea-ice
coverage are common; both inhospitable conditions for most foraminifera
species. The non-existence of foraminifera in ocean sediments are rarely
reported, since the absence of foraminifera leaves little to be analysed in
ecological or geochemical studies. Thus, observations may be biased towards
the warm climate events naturally present within the inter-decadal
variability, when foraminifera are found. We account for this through
specification of measurement error in (6), $\textbf{W}_{\textbf{X}}$. Define
$B_{N}$ as a set of indices that index observations from the Nordic Seas in
$\textbf{Z}_{\textbf{X}}$. We set $\mathbb{E}[\textbf{W}_{\textbf{X}}]_{b}=2$
and $\mathbb{E}[\textbf{W}_{\textbf{X}}]_{b^{*}}=0$ for $b\in B_{N}$ and
$b^{*}\notin B_{N}$. As the bias originates from a systematic reporting error
we believe the measurement errors to be correlated. We partition
$\text{var}[\textbf{W}_{\textbf{X}}]=V_{\textbf{D}}(\textbf{I}+V_{\textbf{B}})V_{\textbf{D}}$
and set $V_{\textbf{D}}$ as the diagonal matrix of the reported standard
deviations in the MARGO dataset, and $V_{\textbf{B}[b,b^{\prime}]}=1\;\forall
b,b^{\prime}\in B_{N}$ where $b\neq b^{\prime}$, and $0$ otherwise. Finally,
we set incidence matrix in (6),
$\textbf{H}_{\textbf{X}}=\textbf{H}_{\textbf{X}}^{T}\otimes\textbf{H}_{\textbf{X}}^{S}$,
where $\textbf{H}_{\textbf{X}}^{T}$ calculates either the annual or summer
means of $\textbf{X}^{*}$ as necessary, and $\textbf{H}_{\textbf{X}}^{S}$
spatially interpolates these averages from simulator’s spatial grid to the
data locations.

####  4.1.3 SIC Observations

Point-wise proxy measurements of SIC are difficult to interpret and
unreliable. More robust are estimates of maximum (winter) sea-ice extent. The
Northern and Southern extents used for $\textbf{Z}_{\textbf{Y}}$ are shown in
Figure 1 by the red dashed lines. The Southern Hemisphere maximum sea-ice
extent is as published in Gersonde et al. (2005). Northern Hemisphere extents
are only available for specific regions (e.g. De Vernal et al. (2005); Crosta
et al. (1998); Xiao et al. (2015)), so we use a simple estimate of the
Northern Hemisphere sea-ice extent provided by our domain experts. We specify
$\textbf{Z}_{\textbf{Y}}$ as a $p$-dimensional vector that records a $1$ at
spatial locations within the extent boundaries, and a $0$ outside. Incidence
matrix, $\textbf{H}_{\textbf{Y}}$ in (14) is a $p\times q$ matrix that
collates, from $\textbf{Y}^{*}(\textbf{X}^{*})$, the February SIC from the
Northern hemisphere and August SIC from the Southern hemisphere. We specify
measurement error, $\textbf{W}_{\textbf{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}[\textbf{W}_{\textbf{Y}}]=K\text{cor}[\textbf{W}_{\textbf{Y}}]K$
where $K$ is a diagonal matrix that represents our marginal uncertainty
(capturing the fact that we are more uncertain of the coverage the nearer we
get to the extent boundary) and $\text{cor}[\textbf{W}_{\textbf{Y}}]$ is
spatially correlated error; descriptions of $K$ and
$\text{cor}[\textbf{W}_{\textbf{Y}}]$ are given in the supplementary material.

###  4.2 Fitting the coexchangeable process model

To fit the coexchangeable process model we first model SST via the
coexchangeable model described in Section 2.2 and the process of SIC given SST
is modelled using the methodology developed in Section 3.

####  4.2.1 Sea-surface temperature

Following notation in Section 2.2, define the ensemble SST’s as
$\mathbb{X}=\\{\textbf{X}_{1},\dots,\textbf{X}_{m}\\}$, with each
$\textbf{X}_{i}$ as a second order exchangeable vector of length $q=12p$. This
leads to the representation in equation (4), for which we require prior
specification of $\text{var}[\mathcal{M}(\textbf{X})]$ and
$\text{var}[\mathcal{R}_{i}(\textbf{X})]$. Coexchangeabilty of
$\textbf{X}_{i}$ and $\textbf{X}^{*}$ leads to (5), for which we require prior
specification of model mismatch terms $\mathbb{E}[\textbf{U}_{\textbf{X}}]$,
$\text{var}[\textbf{U}_{\textbf{X}}]$ and incidence matrix A. We specify
$\text{var}[\mathcal{R}(\textbf{X})]=\alpha^{2}\text{var}[\mathcal{M}(\textbf{X})]$
so that
$\text{var}[\textbf{X}]=(1+\alpha^{2})\text{var}[\mathcal{M}(\textbf{X})]$
and, following Rougier et al. (2013), set
$\mathbb{E}[\textbf{U}_{\textbf{X}}]=\bm{0}$ and $\textbf{A}=\textbf{I}$.

For SST (and most climate fields) 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}[\textbf{X}]$ and $\text{var}[\textbf{X}^{*}]$ (and hence
$\text{var}[\textbf{U}_{\textbf{X}}]$) have a Kronecker structure. For
example,
$\text{var}[\textbf{X}]=\text{var}[\textbf{X}_{T}]\otimes\text{var}[\textbf{X}_{S}]$
where $T$ denotes time and $S$ denotes space. If we similarly partition
$\text{var}[\textbf{U}_{\textbf{X}}]=\text{var}[\textbf{U}_{\textbf{X}_{T}}]\otimes\text{var}[\textbf{U}_{\textbf{X}_{S}}]$
and equate either
$\text{var}[\textbf{X}_{T}]=\text{var}[\textbf{U}_{\textbf{X}_{T}}]$ or
$\text{var}[\textbf{X}_{S}]=\text{var}[\textbf{U}_{\textbf{X}_{S}}]$
computation of our adjusted beliefs is efficient. Note that this assumption is
far weaker than the application in Rougier et al. (2013) where it is simply
assumed that $\kappa^{2}\text{var}[\textbf{X}]=\text{var}[\textbf{U}]$. Here,
we set
$\text{var}[\textbf{X}_{T}]=\text{var}[\textbf{U}_{\textbf{X}_{T}}]=\textbf{J}_{n\times
n}$ as our domain experts believe that discrepancies between the
$\textbf{X}_{i}$ and reality are predominantly spatial and constant in time.

We set $\text{var}[\textbf{X}_{S}]$ as the positive semi-definite matrix that
minimises the distance between $\text{var}[\textbf{X}]$ and the empirical
covariance matrix of $\mathbb{X}$, $\textbf{S}_{\textbf{X}}$, under the
Frobenius norm. This specification follows similar arguments to Rougier et al.
(2013) where $\text{var}[\textbf{X}]=\textbf{S}_{\textbf{X}}$, but preserves
Kronecker structure in $\text{var}[\textbf{X}]$ to allow for efficient
computation. Details of this calculation are in the supplementary material.
Elements of $\text{var}[\textbf{U}_{\textbf{X}_{S}}]$ are defined via the
$C^{4}$-Wendland covariance function such that

$\text{var}[\textbf{U}_{\textbf{X}_{S}}]_{[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}_{+}$
(20)

where $\tau\geq 6$, $c\in(0,\pi]$, $(a)_{+}=\text{max}(0,a)$, and $d(i,j)$ is
the geodesic distance between locations $i$ and $j$. The $C^{4}$-Wendland
covariance function is commonly chosen so as to define a smooth process on the
sphere; (see, for example Astfalck et al., 2019). Parameters are specified as
$\kappa^{2}=1.61$, $c=0.92$, and $\tau=6$; this represents the prior-beliefs
of domain experts as to the behaviour of $\textbf{U}_{\textbf{X}}$.

The 2-stage Bayes linear update follows Rougier et al. (2013) and is given in
the supplementary material. As with Rougier et al. (2013), we assume the first
update $\mathbb{E}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$ is well approximated by
$\mathbb{E}_{\bar{\textbf{X}}}[\textbf{X}^{*}]\approx\bar{\textbf{X}}$ and we
calculate
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]=\left(\frac{\alpha^{2}}{m+\alpha^{2}}\right)\text{var}[\mathcal{M}(\textbf{X})]+\text{var}[\textbf{U}_{\textbf{X}}]$
so
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]\rightarrow\text{var}[\textbf{U}]$
as $m\rightarrow\infty$. Here, we choose $\alpha^{2}=1$ and have $m=13$ and so
do not assume
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]\approx\text{var}[\textbf{U}]$.
From specifications
$\text{var}[\textbf{U}_{\textbf{X}_{T}}]=\text{var}[\textbf{X}_{T}]$ and
$\text{var}[\textbf{X}]=2\text{var}[\mathcal{M}(\textbf{X})]$,
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$ has Kronecker structure so
that
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]=\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{T}]\otimes\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{S}]$,
where $\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{T}]=\textbf{J}_{n\times
n}$ and
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{S}]=\frac{\text{var}[\textbf{X}_{S}]}{2(m+1)}+\text{var}[\textbf{U}_{\textbf{X}_{S}}]$.
The above specifications lead to a prior predictive for $\textbf{X}^{*}$
(prior to observing $\textbf{Z}_{\textbf{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^{\degree}$C). This is problematic due to the data
sparsity in the poles, and so we correct our prior by adding 10 equally spaced
(over longitude) pseudo-observations at 80°N and 80°S, each of
$-1.92^{\degree}$C. Figures 2(a) and 2(b) plot
$\mathbb{E}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$ and marginal variance of
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$, respectively, for January.

The second update is calculated by

$\displaystyle\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$
$\displaystyle=\bar{\textbf{X}}+\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]\textbf{H}_{\textbf{X}}^{\intercal}\text{var}[\textbf{Z}]^{-1}(\textbf{Z}_{\textbf{X}}-\textbf{H}_{\textbf{X}}\bar{\textbf{X}}-\mathbb{E}[\textbf{W}_{\textbf{X}}])\text{,
and}$ (21)
$\displaystyle\text{var}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$
$\displaystyle=\textbf{J}_{n\times
n}\otimes\left(\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{S}]-\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{S}]\left(\textbf{H}_{\textbf{X}}^{S}\right)^{\intercal}\text{var}[\textbf{Z}]^{-1}\textbf{H}_{\textbf{X}}^{S}\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}_{S}]\right),$
(22)

where
$\text{var}[\textbf{Z}]=\textbf{H}_{\textbf{X}}\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]\textbf{H}_{\textbf{X}}^{\intercal}+\text{var}[\textbf{W}_{\textbf{X}}]$.
These updates are shown in Figure 2(c) (for
$\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$) and
Figure 2(d) (for marginal variance of
$\text{var}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$), also
for January. Figure 2(e) plots the contribution of the second update to our
expected beliefs,
$\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}}[\textbf{X}^{*}]-\mathbb{E}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$.
Regions of cooling are apparent westward of large continental masses (i.e. 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.

(a)

(b)

(c)

(d)

(e)

Figure 2: Adjusted beliefs of January SST: (a) adjusted expectation of SST by
$\bar{\textbf{X}}$, $\mathbb{E}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$; (b) the
adjusted variance of SST by $\bar{\textbf{X}}$,
$\text{var}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$; (c) adjusted expectation of
SST by $\bar{\textbf{X}}$ and $\textbf{Z}_{\textbf{X}}$,
$\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$; (d)
the adjusted variance of SST by $\bar{\textbf{X}}$ and
$\textbf{Z}_{\textbf{X}}$,
$\text{var}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$; and
(e) the contribution of the second update to our expected beliefs of SST
$\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]-\mathbb{E}_{\bar{\textbf{X}}}[\textbf{X}^{*}]$.
Temperature is shown in degrees Celsius and plots of variance show marginal
variance.

####  4.2.2 Sea-ice concentration given sea-surface temperature

Following notation in Section 3 we now consider the process of SIC given SST.
Define the jointly-observed ensemble of SST and SIC as
$(\mathbb{X},\mathbb{Y})=\\{(\textbf{X}_{1},\textbf{Y}_{1}),\dots,(\textbf{X}_{m},\textbf{Y}_{m})\\}$
where each $\textbf{Y}_{i}$ is a $12p$-dimensional vector of SIC that we
explicitly model as a function of SST, $\textbf{Y}_{i}(\textbf{X}_{i})$. The
$\textbf{Y}_{i}(\textbf{X}_{i})$ comprise of $12$ $p$-dimensional
conditionally second order exchangeable spatial processes, $Y_{i,j}(X_{i,j})$,
where $i$ indexes the ensemble member and $j$ indexes time. We assume
$\mathbb{E}[Y_{i,j}(X_{i,j})]=\phi(X_{i,j})\beta_{i}$ which leads to the
representation in (8), for which we require specifications of $\phi(X_{i,j})$,
$\mathbb{E}[\beta_{i}]$, $\text{var}[\beta_{i}]$, and
$\text{var}[\mathcal{R}_{i}(Y_{i};X_{i})]$. We assume the $\beta_{i}$ to be
second order exchangeable which leads to the representation in (9) that
requires specification of $\mathbb{E}[\mathcal{M}(\beta)]$,
$\text{var}[\mathcal{M}(\beta)]$ and $\text{var}[\mathcal{R}(\beta)]$.
Together, these exchangeability judgements lead to the representation for the
$\textbf{Y}_{i}$ in (10) (and equivalently (11)).

For modelling SIC given SST (and for most climate processes) 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
$\phi(X_{i,j})=\Psi(X_{i,j})\bm{\Theta}_{i}$ where $\Psi(X_{i,j})$ models the
behaviour of SIC and SST at individual locations using spline bases, and
$\bm{\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 $\phi(X_{i,j})=\Psi(X_{i,j})$, but as the size of $\beta_{i}$ then scales
with the spatial resolution this is not feasible for climate models. Here we
specify $\Psi(X_{i,j})$ with I-spline bases, a basis commonly used for
monotone functions (such as SIC given SST). $\bm{\Theta}_{i}$ is specified
using a principal component design calculated using a projection of the
$\textbf{Y}_{i}$ onto the column space of the $\Psi(\textbf{X}_{i})$,
$\hat{\Theta}_{i}$. Approximating the spatial coefficients using principal
components restricts inference for the $\beta_{i}$ and $\mathcal{M}(\beta)$ to
linear combinations of $\hat{\Theta}_{i}$; we argue this is appropriate for
modelling the mean MME process. We use a more flexible parameterisation in the
model discrepancy below so that inference on $\textbf{Y}^{*}(\textbf{X}^{*})$
given the SIC data is not restricted to linear combinations of the principal
component design. Full specification of $\Psi(X_{i,j})$ and the
$\bm{\Theta}_{i}$ is in the supplementary material. All
$\Phi_{i}=\phi(\textbf{X}_{i})$ are full rank and so the conditions of Theorem
3.1 are met and $\hat{\bm{\beta}}=(\hat{\beta}_{1},\dots,\hat{\beta}_{m})$ is
sufficient for $(\mathbb{X},\mathbb{Y})$ for updating
$\textbf{B}=(\beta_{1},\dots,\beta_{m},\mathcal{M}(\beta))$.

We specify $\text{var}[\mathcal{R}_{i}(Y_{i};X_{i})]$ as a heteroskedastic
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^{\degree}$C and $3^{\degree}$C where sea-ice behaviour is variable.
Similar to Rougier et al. (2013) we set $\text{var}[\beta_{i}]=\lambda$ where
$\lambda$ is a diagonal matrix of the eigenvalues calculated from the
principal component decomposition described above. As before we set
$\text{var}[\mathcal{R}(\beta)]=\alpha^{2}\text{var}[\mathcal{M}(\beta)]$.

Adjusted beliefs $\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{B}]$ and
$\text{var}_{\hat{\bm{\beta}}}[\textbf{B}]$ (and so to
$\mathbb{E}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ and
$\text{var}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$) are calculated as in the
supplementary material. As mentioned in Section 3.1 an advantage of specifying
our model jointly over the $\beta_{i}$ and $\mathcal{M}(\beta)$ is that we
obtain inference for each $\beta_{i}$ (as opposed to only $\mathcal{M}(\beta)$
as in Goldstein and Wooff (1998)). We look at the individual fits in Figure 3
where we plot SIC given SST, parameterised by the $\beta_{i}$, for each MME
member, at four locations. 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 disagrees with the rest of the ensemble (this is
seen over the entire North Atlantic region and is not the result of a local
erroneous fit). Finally, location D shows sea-ice behaviour at the Antarctic
sea-ice edge. Given the SSTs that are observed, the differing physics
predominantly lead to differences in Winter sea-ice.

Figure 3: SIC given SST, parameterised by $\beta_{i}$, for each ensemble
member, at four locations. Shaded regions denote $\pm 2$ standard deviations.

Coexchangeability of $\textbf{Y}^{*}(\textbf{X}^{*})$ and $\beta_{i}$ leads to
(13) for which specification of $A(\textbf{X}^{*})$,
$\mathbb{E}[\textbf{U}_{\textbf{Y}}]$ and
$\text{var}[\textbf{U}_{\textbf{Y}}]$ is required. We set
$A(\textbf{X}^{*})=\phi(\textbf{X}^{*})$ and model discrepancy so that
$\textbf{U}_{\textbf{Y}}=\Psi(\textbf{X}^{*})\textbf{U}_{\Theta}$ (and thus
$\mathbb{E}[\textbf{U}_{\textbf{Y}}]=\Psi(\textbf{X}^{*})\mathbb{E}[\textbf{U}_{\Theta}]$
and
$\text{var}[\textbf{U}_{\textbf{Y}}]=\Psi(\textbf{X}^{*})\text{var}[\textbf{U}_{\Theta}]\Psi(\textbf{X}^{*})^{\intercal}$).
$\textbf{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 $\textbf{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}[\textbf{U}_{\beta}]=0$ and
$\text{var}[\textbf{U}_{\Theta}]=\textbf{I}_{l\times
l}\text{var}[\textbf{U}_{S}]$ where $l$ is the number of spline coefficients
at each location and $\text{var}[\textbf{U}_{S}]$ is calculated by the
$C^{4}$-Wendland function in (20) with parameters $\kappa^{2}=0.3$, $c=4$, and
$\tau=6$. The first update
$\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$ and
$\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$ is calculated
given $\textbf{X}^{*}$ using
$\mathbb{E}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ and
$\text{var}_{\hat{\bm{\beta}}}[\mathcal{M}(\beta)]$ as previously calculated.
The second update
$\mathbb{E}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$
and
$\text{var}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$
are calculated given $\textbf{X}^{*}$ as in (17) and (18). Expectations and
marginal variances of both updates, calculated with
$\textbf{X}^{*}=\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$
are shown in Figure 4. The expected SSTs do not produce adequate sea-ice when
only considering learnt relations from the MME (Figure 4(a)) and marginal
variance of SIC is large in regions where the expected SST is cold enough to
support sea-ice (Figure 4(b)). Updating SIC using the data leads to more
extensive sea-ice cover (Figure 4(c)) and a reduction of variance everywhere
except the sea-ice boundary (Figure 4(d)). Similar behaviour is seen in the
Northern Hemisphere winter.

(a)

(b)

(c)

(d)

Figure 4: Adjusted beliefs of August SIC given expected SST (i.e.
$\textbf{X}^{*}=\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]$):
(a) the adjusted expectation of SIC by $\hat{\bm{\beta}}$,
$\mathbb{E}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$; (b) the
adjusted variance of SIC by $\hat{\bm{\beta}}$,
$\text{var}_{\hat{\bm{\beta}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$; (c) the
adjusted expectation of SIC by $\hat{\bm{\beta}}$ and
$\textbf{Z}_{\textbf{Y}}$,
$\mathbb{E}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$;
and (d) the adjusted variance of SIC by $\hat{\bm{\beta}}$ and
$\textbf{Z}_{\textbf{Y}}$
$\text{var}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\textbf{X}^{*})]$.
Values are of sea-ice concentration and plots of variance show marginal
variance.

##  5 Sampled boundary conditions and their influence on glacial ice sheet
modelling outputs

The remit of this work was to reconstruct, with uncertainty, joint SST and SIC
fields to act as boundary conditions into FAMOUS-Ice an ‘atmosphere only’
global climate model coupled to a ice sheet model. To briefly 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. To
generate a joint sample ($\tilde{\textbf{X}}^{*}$, $\tilde{\textbf{Y}}^{*}$),
we first sample $\tilde{\textbf{X}}^{*}$ and then
$\tilde{\textbf{Y}}^{*}(\tilde{\textbf{X}}^{*})$. As we have only expectations
and variances for these quantities, we appeal to Chebyshev’s inequality to
define a plausible set for
$\\{\textbf{X}^{*},\textbf{Y}^{*}(\textbf{X}^{*})\\}$, and sample this set
uniformly (similar to how ensemble design is considered in history matching
(Salter et al., 2019)). We sample

$\displaystyle\tilde{\textbf{X}}^{*}$
$\displaystyle\sim\mathbb{E}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}]+(\text{var}_{\bar{\textbf{X}},\textbf{Z}_{\textbf{X}}}[\textbf{X}^{*}])^{1/2}Z$
(23) $\displaystyle\tilde{\textbf{Y}}^{*}(\tilde{\textbf{X}}^{*})$
$\displaystyle\sim\mathbb{E}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\tilde{\textbf{X}}^{*})]+(\text{var}_{\hat{\bm{\beta}},\textbf{Z}_{\textbf{Y}}}[\textbf{Y}^{*}(\tilde{\textbf{X}}^{*})])^{1/2}Z$
(24)

where $(A)^{1/2}$ is any square-root decomposition of matrix $A$ and
$Z\sim\mathcal{U}(-3,3)$. An example of a joint sample
$(\tilde{\textbf{Y}}^{*},\tilde{\textbf{X}}^{*})$ for the months of February
and August is given in Figure 5.

(a) SST - February

(b) SST - August

(c) SIC - February

(d) SIC - August

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

We generate an ensemble of 25 runs comprising a single reference run utilising
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 modelling 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
6. The simulator is run for 5000 years, with the ice sheet initialised 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 et al., 2020).

Figure 6: 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 runs with varying SST and SIC fields in
black.

Figure 6 shows that the boundary conditions have a strong influence on the
simulated ice sheet. SST and SIC can affect ice sheets either though 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 PMIP4 model simulations. Indeed the second update causes a
strong Pacific cooling along the coast of North America (Figure 2e) relative
to the expectation 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 modelling 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.

##  6 Discussion

By exploiting natural conditional exchangeability judgements we develop theory
for the coexchangeable process model, as an extension to Rougier et al.
(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, for
example, Sansom et al. (2021)). Our methodology requires fewer assumptions and
less onerous belief specifications than that required by a probabilistic
Bayesian analysis.

To achieve these advancements we developed a Bayes linear analogue to a
hierarchical Bayesian model. By combining exchangeability judgements and using
the reparameterisation 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
shrinkage and 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 modelling, when any
simulation of part of the Earth system, requires other spatio-temporal fields
as boundary conditions. Our application looked at Palaeo ice-sheet modelling,
where our model had a coupled atmosphere and ice-sheet, with the sea-surface
temperature (SST) and sea-ice concentration (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.

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 modelling. We demonstrated its
efficacy by reconstructing last glacial maximum SST and SIC to force an ice-
sheet model. We showed 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
was 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.

##  7 Acknowledgements

All authors were funded by UKRI Future Leaders Fellowship MR/S016961/1.
Climate-ice sheet simulations were undertaken on ARC4, part of the High
Performance Computing Facilities at the University of Leeds.

## References

  * Annan and Hargreaves (2013) Annan, J., Hargreaves, J.C., 2013\.  A new global reconstruction of temperature changes at the Last Glacial Maximum.  Climate of the Past 9, 367–376.
  * Astfalck et al. (2019) Astfalck, L., Cripps, E., Gosling, J., Milne, I., 2019\.  Emulation of vessel motion simulators for computationally efficient uncertainty quantification.  Ocean Engineering 172, 726–736.
  * Barthelmé and Chopin (2014) Barthelmé, S., Chopin, N., 2014\.  Expectation propagation for likelihood-free inference.  Journal of the American Statistical Association 109, 315–333.
  * Benz et al. (2016) Benz, V., Esper, O., Gersonde, R., Lamy, F., Tiedemann, R., 2016.  Last Glacial Maximum sea surface temperature and sea-ice extent in the Pacific sector of the Southern Ocean.  Quaternary Science Reviews 146, 216–237.
  * Bissiri et al. (2016) Bissiri, P.G., Holmes, C.C., Walker, S.G., 2016.  A general framework for updating belief distributions.  Journal of the Royal Statistical Society. Series B, Statistical methodology 78, 1103.
  * Chandler (2013) Chandler, R.E., 2013.  Exploiting strength, discounting weakness: combining information from multiple climate simulators.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120388\.
  * CLIMAP Project (1981) CLIMAP Project, 1981.  Seasonal reconstructions of the Earth’s surface at the last glacial maximum.  Geological Society of America.
  * Cressie and Johannesson (2008) Cressie, N., Johannesson, G., 2008\.  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70, 209–226.
  * Crosta et al. (1998) Crosta, X., Pichon, J.J., Burckle, L., 1998.  Application of modern analog technique to marine Antarctic diatoms: Reconstruction of maximum sea-ice extent at the last glacial maximum.  Paleoceanography 13, 284–297.
  * De Finetti (1975) De Finetti, B., 1975.  Theory of probability: A critical introductory treatment. volume 1.  John Wiley & Sons.
  * De Vernal et al. (2005) De Vernal, A., Eynaud, F., Henry, M., Hillaire-Marcel, C., Londeix, L., Mangin, S., Matthießen, J., Marret, F., Radi, T., Rochon, A., et al., 2005\.  Reconstruction of sea-surface conditions at middle to high latitudes of the Northern Hemisphere during the Last Glacial Maximum (LGM) based on dinoflagellate cyst assemblages.  Quaternary Science Reviews 24, 897–924.
  * Eyring et al. (2019) Eyring, V., Cox, P.M., Flato, G.M., Gleckler, P.J., Abramowitz, G., Caldwell, P., Collins, W.D., Gier, B.K., Hall, A.D., Hoffman, F.M., et al., 2019\.  Taking climate model evaluation to the next level.  Nature Climate Change 9, 102–110.
  * Gersonde et al. (2005) Gersonde, R., Crosta, X., Abelmann, A., Armand, L., 2005\.  Sea-surface temperature and sea ice distribution of the Southern Ocean at the EPILOG Last Glacial Maximum—a circum-Antarctic view based on siliceous microfossil records.  Quaternary Science Reviews 24, 869–896.
  * Goldstein (1986) Goldstein, M., 1986.  Exchangeable belief structures.  Journal of the American Statistical Association 81, 971–976.
  * Goldstein and Wooff (2007) Goldstein, M., Wooff, D., 2007\.  Bayes linear statistics: Theory and methods. volume 716\.  John Wiley & Sons.
  * Goldstein and Wooff (1998) Goldstein, M., Wooff, D.A., 1998\.  Adjusting exchangeable beliefs.  Biometrika 85, 39–54.
  * Gregoire et al. (2016) Gregoire, L.J., Otto-Bliesner, B., Valdes, P.J., Ivanovic, R., 2016\.  Abrupt Bølling warming and ice saddle collapse contributions to the Meltwater Pulse 1a rapid sea level rise.  Geophysical research letters 43, 9130–9137.
  * Gregoire et al. (2012) Gregoire, L.J., Payne, A.J., Valdes, P.J., 2012.  Deglacial rapid sea level rises caused by ice-sheet saddle collapses.  Nature 487, 219\.
  * Gregory et al. (2020) Gregory, J.M., George, S.E., Smith, R.S., 2020.  Large and irreversible future decline of the greenland ice sheet.  The Cryosphere 14, 4299–4322.
  * Harrison et al. (2015) Harrison, S.P., Bartlein, P., Izumi, K., Li, G., Annan, J., Hargreaves, J., Braconnot, P., Kageyama, M., 2015\.  Evaluation of CMIP5 palaeo-simulations to improve climate projections.  Nature Climate Change 5, 735–743.
  * Hodges (1998) Hodges, J.S., 1998.  Some algebra and geometry for hierarchical models, applied to diagnostics.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 60, 497–536.
  * Hoff (2009) Hoff, P.D., 2009.  A first course in Bayesian statistical methods. volume 580.  Springer.
  * Ivanovic et al. (2016) Ivanovic, R., Gregoire, L., Kageyama, M., Roche, D., Valdes, P., Burke, A., Drummond, R., Peltier, W.R., Tarasov, L., 2016.  Transient climate simulations of the deglaciation 21-9 thousand years before present; PMIP4 core experiment design and boundary conditions.  Geoscientific Model Development 9, 2563–2587.
  * Jesus and Chandler (2011) Jesus, J., Chandler, R.E., 2011\.  Estimating functions and the generalized method of moments.  Interface focus 1, 871–885.
  * Kageyama et al. (2017) Kageyama, M., Albani, S., Braconnot, P., Harrison, S.P., Hopcroft, P.O., Ivanovic, R.F., Lambert, F., Marti, O., Peltier, W.R., Peterschmitt, J.Y., et al., 2017.  The PMIP4 contribution to CMIP6–Part 4: Scientific objectives and experimental design of the PMIP4-CMIP6 Last Glacial Maximum experiments and PMIP4 sensitivity experiments.  Geoscientific Model Development 10, 4035–4055.
  * Kageyama et al. (2021) Kageyama, M., Harrison, S.P., Kapsch, M.L., Lofverstrom, M., Lora, J.M., Mikolajewicz, U., Sherriff-Tadano, S., Vadsaria, T., Abe-Ouchi, A., Bouttes, N., et al., 2021\.  The PMIP4 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3 simulations.  Climate of the Past 17, 1065–1089.
  * Kucera et al. (2005) Kucera, M., Rosell-Melé, A., Schneider, R., Waelbroeck, C., Weinelt, M., 2005.  Multiproxy approach for the reconstruction of the glacial ocean surface (MARGO).  Quaternary Science Reviews 24, 813–819.
  * Liu and Guillas (2017) Liu, X., Guillas, S., 2017.  Dimension reduction for Gaussian process emulation: an application to the influence of bathymetry on tsunami heights.  SIAM/ASA Journal on Uncertainty Quantification 5, 787–812.
  * Masson and Knutti (2011) Masson, D., Knutti, R., 2011\.  Climate model genealogy.  Geophysical Research Letters 38.
  * Paul et al. (2020) Paul, A., Mulitza, S., Stein, R., Werner, M., 2020\.  A global climatology of the ocean surface during the Last Glacial Maximum mapped on a regular grid (GLOMAP).  Climate of the Past Discussions , 1–26.
  * Rougier et al. (2013) Rougier, J., Goldstein, M., House, L., 2013.  Second-order exchangeability analysis for multimodel ensembles.  Journal of the American Statistical Association 108, 852–863.
  * Salter et al. (2018) Salter, J.M., Williamson, D.B., Gregoire, L.J., Edwards, T.L., 2018\.  Quantifying spatio-temporal boundary condition uncertainty for the North American deglaciation.  arXiv preprint arXiv:1808.09322 .
  * Salter et al. (2019) Salter, J.M., Williamson, D.B., Scinocca, J., Kharin, V., 2019\.  Uncertainty quantification for computer models with spatial output using calibration-optimal bases.  Journal of the American Statistical Association 114, 1800–1814.
  * Sansom et al. (2021) Sansom, P.G., Stephenson, D.B., Bracegirdle, T.J., 2021.  On constraining projections of future climate using observations and simulations from multiple climate models.  Journal of the American Statistical Association , 1–12.
  * Sarnthein et al. (2003) Sarnthein, M., Gersonde, R., Niebler, S., Pflaumann, U., Spielhagen, R., Thiede, J., Wefer, G., Weinelt, M., 2003\.  Overview of glacial Atlantic Ocean mapping (GLAMAP 2000).  Paleoceanography 18.
  * Schmidt et al. (2014) Schmidt, G., Annan, J., Bartlein, P., Cook, B., Guilyardi, É., Hargreaves, J., Harrison, S., Kageyama, M., Legrande, A.N., Konecky, B., et al., 2014\.  Using palaeo-climate comparisons to constrain future projections in CMIP5.  Climate of the Past 10, 221–250.
  * Sha et al. (2019) Sha, Z., Rougier, J., Schumacher, M., Bamber, J., 2019\.  Bayesian model–data synthesis with an application to global glacio-isostatic adjustment.  Environmetrics 30, e2530.
  * Smith et al. (2021) Smith, R.S., George, S., Gregory, J.M., 2021.  FAMOUS version xotzt (FAMOUS-ice): a general circulation model (GCM) capable of energy-and water-conserving coupling to an ice sheet model.  Geoscientific Model Development 14, 5769–5787.
  * Tarasov et al. (2012) Tarasov, L., Dyke, A.S., Neal, R.M., Peltier, W.R., 2012\.  A data-calibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling.  Earth and Planetary Science Letters 315, 30–40.
  * Tierney et al. (2020) Tierney, J.E., Zhu, J., King, J., Malevich, S.B., Hakim, G.J., Poulsen, C.J., 2020\.  Glacial cooling and climate sensitivity revisited.  Nature 584, 569–573.
  * Whitaker and Hamill (2002) Whitaker, J.S., Hamill, T.M., 2002\.  Ensemble data assimilation without perturbed observations.  Monthly weather review 130, 1913–1924.
  * Williamson and Sansom (2019) Williamson, D.B., Sansom, P.G., 2019\.  How are emergent constraints quantifying uncertainty and what do they leave behind?  Bulletin of the American Meteorological Society 100, 2571 – 2588.
  * Xiao et al. (2015) Xiao, X., Fahl, K., Müller, J., Stein, R., 2015\.  Sea-ice distribution in the modern Arctic ocean: Biomarker records from trans-Arctic ocean surface sediments.  Geochimica et Cosmochimica Acta 155, 16–29.
  * Zhang et al. (2020) Zhang, B., Cressie, N., et al., 2020\.  Bayesian inference of spatio-temporal changes of Arctic sea ice.  Bayesian Analysis 15, 605–631.