Paleoclimatic and paleoenvironmental reconstructions are fundamentally uncertain because no proxy is a direct record of a single environmental variable of interest; all proxies are indirect and sensitive to multiple forcing factors. One productive approach to reducing proxy uncertainty is the integration of information from multiple proxy systems with complementary, overlapping sensitivity. Mostly, such analyses are conducted in an ad hoc fashion, either through qualitative comparison to assess the similarity of single-proxy reconstructions or through step-wise quantitative interpretations where one proxy is used to constrain a variable relevant to the interpretation of a second proxy. Here we propose the integration of multiple proxies via the joint inversion of proxy system and paleoenvironmental time series models in a Bayesian hierarchical framework. The “Joint Proxy Inversion” (JPI) method provides a statistically robust approach to producing self-consistent interpretations of multi-proxy datasets, allowing full and simultaneous assessment of all proxy and model uncertainties to obtain quantitative estimates of past environmental conditions. Other benefits of the method include the ability to use independent information on climate and environmental systems to inform the interpretation of proxy data, to fully leverage information from unevenly and differently sampled proxy records, and to obtain refined estimates of proxy model parameters that are conditioned on paleo-archive data. Application of JPI to the marine

Paleoenvironmental reconstructions, including reconstructions of past
climate, provide a powerful tool to document the sensitivity of Earth
systems to forcing, characterize the range of natural responses associated
with different modes of global change, and identify key mechanisms governing
these responses. Throughout the vast majority of the planet's history,
however, estimates of environmental conditions can only be obtained through
proxy reconstructions. The word proxy is derived from the Latin word

The simplest proxy reconstructions typically focus on a single environmental variable of interest. Experimental or natural calibration datasets are used to calibrate mathematical relationships between the environmental variable and proxy measure, and these relationships are inverted to obtain quantitative estimates of that variable. Residual variance in the calibration is treated as noise. In reality, however, no proxy exists that is sensitive only to a single paleoenvironmentally relevant variable, and a large part of the proxy system noise reflects the uncharacterized influence of other environmental and post-depositional variables. Fossil leaf assemblages, for example, exhibit variability that can be associated with mean annual air temperature but also may be influenced by many other environmental variables and evolutionary history (Royer et al., 2005; Greenwood et al., 2004). The saturation state of alkenones produced by marine phytoplankton is a sensitive recorder of water temperature, but characteristics of alkenones preserved in marine sediments are also strongly affected by physiological factors, seasonality of production, and selective degradation (Conte et al., 1998, 2006). Even recently emerging clumped isotope techniques, which are in theory a direct recorder of the temperature of carbonate mineral formation, can be affected by factors such as growth rate, carbonate system disequilibrium, and poorly constrained, potentially variable offsets between the environment of carbonate formation and more commonly targeted atmospheric temperature conditions (Passey et al., 2010; Affek et al., 2014; Saenger et al., 2012).

Failure to recognize and consider the sensitivity of proxies to multiple environmental factors leads to two important problems in traditional proxy interpretations. First, considering only a single environmental variable in our interpretations maximizes the uncertainty in our reconstructions. Uncertainty could be reduced if the influence of other variables is described and constrained. Second, unacknowledged sensitivity to multiple variables creates potential for biased proxy interpretations if variation in these variables is non-random across the reconstruction.

A productive approach to addressing these issues is the use of proxy system models in the interpretation of proxy data (Evans et al., 2013). These models represent an attempt to mathematically describe the complex of environmental, physical, and biological factors that control how environmental signals are sampled, recorded, and preserved in proxy measurements. Recent reviews and perspectives are available discussing the concepts underlying proxy system models and different ways that they have been applied to proxy interpretation, ranging from substitution for empirical calibrations in inverse estimation of environmental signals to formal integration within climate model data assimilation schemes (Evans et al., 2013; Dee et al., 2016). A growing number of proxy system models and modeling systems are being developed (e.g., Tolwinski-Ward et al., 2011; Stoll et al., 2012; Dee et al., 2015), and useful models span a range of complexity from empirically constrained regressions to mechanistic, theory-based formulations. Key to any such model is accurate representation of uncertainty in each model component, which allows even relatively simple, potentially incomplete models to be used to obtain reconstructions with quantifiable uncertainty bounds.

Reducing the uncertainty of quantitative paleoenvironmental reconstructions, however, further requires adding constraints to proxy interpretations. In situations where two or more proxies share sensitivity to common or complementary environmental variables, it stands to reason that the information provided by each can be used to refine interpretation of the multi-proxy suite. In practice, a variety of approaches have been used. Commonly, multi-proxy integration has been qualitative and focused on confirmation: trends reconstructed using one proxy system are cross-checked against a second, providing increased confidence in the reconstruction where the patterns match and prompting further investigation where they do not (e.g., Grauel et al., 2013; Keating-Bitonti et al., 2011; Zachos et al., 2006). In other cases, proxies have been combined quantitatively, but usually in a stepwise fashion: one proxy system is used to reconstruct an environmental variable to which it is sensitive, and those reconstructed values are then used to constrain the interpretation of a second proxy (e.g., Fricke et al., 1998; Lear et al., 2000). Although it provides a simple strategy to combining complementary proxy information, this approach does not fully leverage overlapping information that may be contained in multiple systems that respond to common forcing, is not conducive to robust quantification of uncertainty, and requires that both proxies sample coeval paleoenvironmental conditions.

Here we propose a general approach to proxy interpretation that leverages
the benefits of proxy models and provides a robust statistical basis for
multi-proxy integration. The method, which we call Joint Proxy Inversion (JPI), couples proxy models with simple environmental time series models
representing paleoenvironmental target variables in a Bayesian hierarchical
modeling framework (Fig. 1). The hierarchical model is then inverted using
Markov Chain Monte Carlo methods (Geman and Geman, 1984) to obtain posterior parameter estimates and paleoenvironmental time series that are conditioned simultaneously on all proxy and calibration data. Similar approaches have been applied to conduct large-scale meta-analyses (Tingley and Huybers, 2010; Li et al., 2010; Tingley et al., 2012; Garreta et al., 2010) but have not found widespread use in quantitative proxy interpretation. We begin by describing an implementation of JPI for the widely used foraminiferal

Implementation of JPI for the coupled

Proxy and proxy model calibration datasets were compiled from published work
(Fig. 1). Estimates from fluid inclusions, calcite veins, large foraminifera, and echinoderm fossils (Dickson, 2002; Coggon et al., 2010; Lowenstein et al., 2001; Evans et al., 2018; Horita et al., 2002) were combined with information on modern seawater

Foraminiferal

Calibration datasets were compiled to constrain the

For

The age of each pre-modern datum was taken from the primary source. Age uncertainties, where known, can be incorporated in the JPI analysis framework by treating ages as random variables rather than as fixed values and/or including proxy model components representing processes governing the time integration of observations. For simplicity, we do not include such a treatment here. In the discussion we note examples where including age uncertainty would produce a more robust analysis.

The proxy system models comprise the “data model” layer of the hierarchical
model, representing how environmental signals are embedded in the
paleo-proxy and proxy calibration data. The models used here are comprised
of simple transfer functions relating proxy data to contemporaneous
environmental variables and as such can be considered “sensor models” in
the terminology of Evans et al. (2013), with aspects of proxy signal integration and sampling treated in the “archive” and “observation” models of those authors being swept into the error terms of our data model Eqs. (1)–(3). The simplest model is that for seawater

We model foraminiferal

Foraminiferal calibration and proxy

Although not treated as such in most reconstructions, paleoenvironmental conditions are autocorrelated in time, meaning that each proxy observation provides information about conditions not just at a single point in time but across a segment of time. To reflect this, we model paleoenvironmental variables as time series using a correlated random walk model. This parameterization is desirable in that it is minimally prescriptive (i.e., no preferred state or pattern of change is proscribed) but allows incorporation of constraints on (and extraction of inference about) two basic characteristics of the paleoenvironmental system – namely its rate and directedness of change. The environmental models represent the “process model” layer of the Bayesian hierarchical model.

The correlated random walk for variable

For seawater

We select the bounds, base resolution, and prior distributions for the
bottom water temperature and

The model structure described above was coded in the BUGS (Bayesian inference Using Gibbs Sampling) language (Lunn et al., 2012), and Markov Chain Monte Carlo was used to generate samples from the posterior distribution of all model parameters conditioned on the proxy and calibration datasets. The analysis was implemented in R version 3.5.1 (R Core Team, 2019) using the rjags (Plummer, 2018) and R2jags (Su and Yajima, 2015) packages. Three to nine chains were run in parallel. Convergence was assessed visually via trace plots and with reference to the Gelman and Rubin convergence factor (Rhat; Gelman and Rubin, 1992) and effective sample sizes reported by rjags.

For the site 806 analysis, nine chains were run to a length of

For the Pleistocene data we conducted three different analyses, the first
two inverting data from each site independently and the third inverting both
records together. For the joint inversion of both records, we treated each
paleoenvironmental time series as independent, i.e., no correlation structure
was imposed on or fit to the conditions simulated at the two sites, and the
model consists of four time series process models (one each for BWT and

Run times for all analyses can be substantially reduced by adopting a smaller number of time steps (e.g., only the base series) and using interpolation to estimate environmental parameter values at the proxy observation time points. Results from experiments using this approach (not shown) were not detectably different from those shown here.

The paleoenvironmental reconstructions obtained by applying JPI to the site 806 data are similar, to 1st order, to the reconstructions from Lear et al. (2015; hereafter L15) on which our analysis was modeled (Figs. 2 and 3). Our estimates of seawater

Reconstructed seawater

Reconstructed bottom water temperature

JPI paleoenvironmental time series for the single- and multi-site analyses
of the Pleistocene data were nearly identical, with slightly broader
credible intervals for both parameters (BWT and

Reconstructed bottom water temperature

We will now examine several characteristics of the paleoenvironmental time
series obtained in the JPI posterior sample and contrast them with
reconstructions obtained through traditional proxy interpretation methods.
One visually striking difference between the JPI and L15 reconstructions is
the higher BWT and

JPI, in contrast, explicitly considers temporal autocorrelation of the
underlying environmental variables, treating each proxy observation as a
sample arising from one or more underlying, autocorrelated environmental
time series. The properties of the time series themselves, rather than being
assumed, are estimated using the proxy models and the data, meaning that the
smoothed reconstruction reflects the information content of the data. For
very certain proxy models or densely distributed data that record
high-frequency variability, the reconstructed time series will express
short-term changes in the environment. In contrast, reconstructions based on
uncertain models or sparsely sampled data will tend toward greater smoothing
and reflect the longer-term evolution of the mean state of the system. This
is nicely illustrated by comparison of JPI

Another advantage of embedding time series models in JPI is that it offers
an explicit framework for integration of differently sampled proxy records.
In most of the studies reviewed here foraminiferal

A final outgrowth of the integration of proxy system and paleoenvironmental
time series models via JPI is that the method provides quantitative
uncertainty bounds that are linked to and reflect the stratigraphic
distribution and density of proxy information. Because environmental
parameters are modeled as continuous time series, estimates of central
tendency and dispersion (e.g., credible intervals) are obtained throughout
the reconstruction period. For time steps in which no observational data are
available, the dispersion of posterior estimates increases consistent with
the properties of the time series model (e.g., between

In addition to estimating the paleoenvironmental record, JPI provides
posterior estimates of parameters in the underlying paleoenvironmental time
series models and proxy (calibration) models, and these themselves can be
informative. Bayesian inversion has previously been used to estimate proxy
model parameter values in situations where these are poorly constrained
(Tolwinski-Ward et al., 2013), and the joint inversion of proxy and environmental time series models performed in JPI can similarly be used to provide constraints on parameter values for all model components (e.g., Fig. S4). Because the proxy system models used here are simple, and the calibration data themselves are used to generate prior estimates on model parameters, the posterior estimates are generally quite similar to the priors (Fig. 5). The only notable exception is

Prior (black) and posterior (red) distributions for

These refinements reflect a combination of the constraints offered by the
calibration and down-core proxy data. Although at first consideration the
relevance of the latter to calibrating proxy model parameters might not be
apparent, in fact the proxy model must not only be consistent with the calibration data but also explain the observed proxy data given the “true”
environmental conditions. As a result, for a given set of proxy data and
environmental time series model properties only a subset of proxy model
parameter values will be plausible. Consider, for example, the proxy model
precision parameter. In our model construction, this value explains the
“noise” both within the model calibration dataset and the proxy record,
each of which can arise from a similar ensemble of factors (e.g., temporal
variation in the environment at timescales below the time series model time
step, biological or random variation in the environment–proxy relationship).
Our analysis suggests that before the mid-Pleistocene transition, the proxy
model variance implied by the full JPI inversion is similar to that
estimated from the calibration data alone (solid curves in Fig. 5d and h),
with slightly higher

Because the JPI analysis involves sampling of all model parameters
simultaneously, it also can identify and account for correlation among
parameters. The proxy model parameter estimates for site 806 provide a clear
example (Fig. 6). The posterior distributions show strong correlation
between the seawater

Bivariate density plots of the posterior distributions for

JPI also provides posterior estimates on the environmental time series model
parameters, and these distributions can provide information complementary to
the reconstructed time series themselves. Comparing prior and posterior
estimates at all three study sites (Fig. 7), the analysis provides strong
posterior constraints on the error autocorrelation (i.e., directedness of
change). Posterior estimates of the error variance (i.e., magnitude of change
between time steps) for

Prior (black) and posterior (red) parameter distributions for bottom water temperature (BWT, solid) and seawater

In contrast, the error autocorrelation term, which reflects the directedness
of environmental change across model time steps, shows substantial variation
among the data sets (Fig. 7, left column). The highest posterior values
(mean values of 0.77 and 0.92 for BWT and

In this final section, we explore additional examples of how JPI results might be used to support inference or hypothesis testing in paleoenvironmental reconstruction. The multivariate posterior samples produced by JPI provide a sound basis for testing hypotheses of change within or between proxy records. Consider the case where we want to assess the magnitude of change in site 806 bottom water temperature relative to the modern (core top) value. Unlike the raw proxy data or traditional interpretations thereof, the JPI samples provide distributions for the environmental variables that support testing at any point in time represented in the paleoenvironmental time series. Other interpolation or smoothing methods can and have been used to conduct such tests, for example of change in global temperature relative to modern (Marcott et al., 2013), but an advantage of JPI, again, is that correlation among model parameters and temporal autocorrelation are included and optimized in the analysis, reducing the need to independently and subjectively specify these.

The effect of parameter correlation can be seen in comparing change relative
to modern within individual posterior samples (within sample) versus change
between each posterior sample and the 0 Ma median value (between sample;
Fig. 8a); the latter being equivalent to a traditional test for non-zero
difference that assumes independence. At short time lags (less than

Evaluating changes within and between environmental
reconstructions using JPI output.

Another example involves cross-site comparison. Here, we similarly ask
whether seawater

Finally, because JPI results provide integrated, self-consistent estimates
of multiple environmental variables, it can be used to identify and
characterize multivariate modes of environmental change in Earth's past.
Results from the site 806 analysis, for example, demonstrate non-linear
coupling between changes in BWT and

Bivariate density plot of posterior values from the site 806 environmental time series models (base 50 kyr time steps only). All values are plotted as change relative to 18 Ma within an individual posterior sample. Dots show the median values from the posterior time series.

Traditional approaches to proxy interpretation suffer from broad and poorly characterized uncertainty and potential biases related to the sensitivity of proxies to multiple environmental factors (Sweeney et al., 2018). Proxy system modeling and multi-proxy reconstruction provide partial solutions to these issues, but a robust accessible framework for integrating these two approaches in the development of paleoenvironmental reconstructions is also needed. We suggest that Bayesian hierarchical models that leverage simple time series representations of paleoenvironmental conditions offer such a framework. This approach is broadly generalizable to any set of proxies for which appropriate forward models can be written. It confers many of the advantages of more complex data assimilation methods that leverage Earth system models (Evans et al., 2013), while remaining independent of the assumptions embedded in these models and flexible enough to be applied over a wide range of systems and timescales. As with any statistically based analysis, JPI results are model-dependent: they provide a basis for interpreting data in the context of a specific model and its assumptions, and this dependence should be acknowledged and considered in the presentation and interpretation of results.

Our illustration of the method based on the coupled

All data and code used to conduct the analyses and create figures reported in this paper are archived online (Bowen, 2019) and available at

The supplement related to this article is available online at:

GJB conceived, designed, and conducted the analyses with support from BFF, AS, and GJR. CHL provided access to data and advice on application of the

The authors declare that they have no conflict of interest.

This research has been supported by the National Science Foundation, Division of Earth Sciences (grant no. 1502786).

This paper was edited by Helen McGregor and reviewed by two anonymous referees.