PaleoSTeHM v1.0: a modern, scalable spatiotemporal hierarchical modeling framework for paleo-environmental data

Geological records of past environmental change provide crucial insights into long-term climate variability, trends, non-stationarity, and nonlinear feedback mechanisms. However, reconstructing spatiotemporal fields from these records is statistically challenging due to their sparse, indirect, and noisy nature. Here, we present PaleoSTeHM, a scalable and modern framework for spatiotemporal hierarchical modeling of paleo-environmental data. This framework enables the implementation of flexible statistical models that rigorously quantify spatial and temporal variability from geological data while clearly distinguishing measurement and inferential uncertainty from process variability. We illustrate its application by reconstructing temporal and spatiotemporal paleo-sea-level changes across multiple locations. Using various modeling and analysis choices, PaleoSTeHM demonstrates the impact of different methods on inference results and computational efficiency. Our results highlight the critical role of model selection in addressing specific paleo-environmental questions, showcasing the PaleoSTeHM framework's potential to enhance the robustness and transparency of paleo-environmental reconstructions.

1 Introduction

As humans push the planet's climate and biosphere increasingly far outside the range of our species' experience, the environmental reconstructions derived from the geological record provide critical out-of-sample information to test the physical models used to project future environmental change. However, as environmental records, geological data are sparse, often noisy, and indirect (PAGES2k Consortium, 2017; Shennan, 2015). Reconstructing paleo-environmental fields is thus a critical and challenging statistical task (Tingley et al., 2012).

From a modeling perspective, spatiotemporal hierarchical statistical models provide a natural, conceptually straightforward framework for reconstructing paleo-environmental signals (Ashe et al., 2019; Cressie and Wikle, 2015; Tingley et al., 2012). Hierarchical statistical models, often employed within a Bayesian framework, decompose the various sources of random variation contributing to individual observations into distinct levels, thereby providing a clear articulation of the assumptions underlying the statistical analysis. They have been increasingly used to model paleo-climate fields from geological proxies, which are naturally occurring physical characteristics or chemical markers that can be used to reconstruct past climate and environmental conditions, such as temperature and precipitation, from sources like tree rings and corals (Walter et al., 2022; PAGES2k Consortium, 2017). These applications have proven crucial in assessing the robustness of scientific knowledge of past climate and placing changes in the modern, instrumentally observed period in the context of longer-term variability. For example, they have shown an increasing influence of ice melt and thermal expansion on the global mean sea level (GMSL) since 1860 CE (Walker et al., 2021), that GMSL rise over the 20th century was faster than during any century in at least 3000 years (Kemp et al., 2018; Kopp et al., 2016), and that several early 21st-century Arctic summers exhibited warmth unprecedented in at least 600 years (Tingley and Huybers, 2013).

Hierarchical models are in high demand within the paleo-environmental research community. For example, in the past few years, numerous papers have used temporal or spatiotemporal hierarchical models with Gaussian process (GP) priors to interpret paleo-sea-level proxies (e.g. Tan et al., 2023; Khan et al., 2022; Vacchi et al., 2021). To meet the demand of the paleo-environmental research community, this paper describes PaleoSTeHM v1.0, which is designed to support the flexible and high-performance implementation of spatiotemporal hierarchical modeling for paleo-environmental data. PaleoSTeHM (https://github.com/radical-collaboration/PaleoSTeHM, last access: 1 April 2025) is a framework built in the spirit of open science and utilizes modern machine learning architecture (e.g. Pollack et al., 2024). It is designed so users can select not only various modeling choices, such as change-point models for temporal analysis or GP for spatiotemporal analysis, but also analysis choices, including fully Bayesian, empirical Bayesian, and variational Bayesian analysis (more details in Sect. 2), to investigate different research questions, with different types of data and spatiotemporal scales (e.g., local to global, years to millennia) considered. In this paper, some key terms and phrases are defined in Table 1.

Table 1 Definitions of relevant terms in this study. This paper employs terminology based on Ashe et al. (2019).

2 Hierarchical statistical modeling

Hierarchical modeling is a statistical approach that separates multiple sources of variability contributing to individual observations into distinct levels, enabling a clear understanding and quantification of uncertainties. This section briefly describes basic theory of hierarchical modeling in the paleo-environment, using paleo-sea level as an illustrative example. For more systematic introductions to hierarchical statistical modeling of paleo-sea level and paleo-climate, readers can refer to Ashe et al. (2019) and Tingley et al. (2012).

Bayesian statistics denotes a statistical theory that uses Bayes' theorem to update probabilities conditioned on data and prior knowledge. Based on Bayes' theorem (Laplace, 1810), the conditional probability of the observed data (y) can be derived from the conditional probability of unknown parameter(s) or process(es) (θ):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }p\left(\mathit{\theta }\right|y)={\displaystyle \frac{p\left(y\right|\mathit{\theta }\left)p\right(\mathit{\theta })}{p\left(y\right)}},\end{array}$$

where p denotes “probability” and | represents “given”. The likelihood function, p(y|θ), represents the probability of observing the data y given the parameter(s) or process(es) θ of the model. The prior distribution, p(θ), captures a priori beliefs about the unknown parameter(s) or process(es) before any data are observed. The term p(y), known as the marginal likelihood (or evidence), is the probability of the observed data averaged over all possible parameters or processes:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }p\left(y\right)=\int p(y\mid \mathit{\theta })p\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }.\end{array}$$

Given the observations, the posterior distribution, p(θ|y), reflects the updated beliefs about the parameter(s) or process(es) after considering the data. Since the marginal likelihood p(y) is often intractable and remains constant for a given dataset, we use the simplified form of Bayes' theorem, where the posterior distribution is proportional to the product of the likelihood and the prior:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }p\left(\mathit{\theta }\right|y)\propto p(y\left|\mathit{\theta }\right)p\left(\mathit{\theta }\right),\end{array}$$

where ∝ indicates “is proportional to”.

A basic hierarchical statistical model distinguishes the change in observations from both its inherent variability and the observational noise. These models achieve probabilistic uncertainty estimation for time series and/or spatial fields by treating observed data as conditional on a latent (unobserved) process and unknown parameters, enabling separate quantification of uncertainties at each level through the application of Bayesian conditional probabilities. Each level of the model quantifies uncertainties independently, necessitating careful evaluation of their respective sources. Generally, three levels are defined: the data level, the process level, and the parameter level.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }p(f,{\mathit{\theta }}_s,{\mathit{\theta }}_d\mid y)\propto \underset{\text{data model}}{\underbrace{p(y\mid f,{\mathit{\theta }}_d)}}\cdot \underset{\text{process model}}{\underbrace{p(f\mid {\mathit{\theta }}_s)}}\cdot \underset{\text{parameter model}}{\underbrace{p({\mathit{\theta }}_d,{\mathit{\theta }}_s)}}\end{array}$$

Taking paleo-relative sea-level (RSL) change as an example, the data-level model defines the relationship between the latent (unobserved) RSL process (f) and the observed RSL data (instrumental and/or proxy), y, while accounting for measurement, inferential (e.g., arising from converting a proxy's elevation to a distribution of RSL), and dating uncertainties (often inherited from geochronology techniques; Reimer et al., 2020; Wright et al., 2017). This level represents the probability distribution of observing a particular sea-level height at a given age, conditioned on the underlying latent process and the associated uncertainties, encapsulated by the data-level parameters, θ_d.

The process level distinguishes the underlying phenomenon of interest and its inherent variability from the noisy observation captured at the data level. This model integrates scientific understanding and associated uncertainties into the estimation of the true RSL process conditioned on model parameters, θ_s. These parameters may represent unobserved physical model parameters (e.g., Earth's rheology in a glacial isostatic adjustment (GIA) model), statistical model parameters (such as the rate of change in a sea-level model), or hyperparameters (parameters of a prior distribution, such as length scale and variance in a GP model). At the foundational level, the parameter model specifies the prior distribution for all unknown parameters, effectively capturing the essential characteristics of both the data and process levels through the unobserved parameters.

In addition to constructing models at the data, process, and parameter levels, often referred to as modeling choices (Ashe et al., 2019), it is essential to choose an appropriate analysis choice for a specific model. This involves decisions regarding the implementation of a model structure, such as deterministic methods or probabilistic methods like Bayesian analysis. Deterministic methods, such as least-squares analysis (Wilks, 1938) and likelihood maximization (Aitken, 1936), rely on fixed relationships between states and events without incorporating randomness into the modeling process. In contrast, probabilistic methods, like Bayesian analysis (Hastings, 1970), account for uncertainty explicitly by representing model parameters and outputs as probability distributions, enabling flexible and robust uncertainty quantification. Analysis choices are also integral to addressing how measurement uncertainties, particularly those arising from geochronological techniques (i.e., input uncertainty), are incorporated and managed within the model. This ensures that the uncertainty is properly quantified and reflected in the final analysis outputs (Ashe et al., 2019). The selection of modeling and analytical choices should consider the problem's complexity, data size and resolution, computational resources, and prior knowledge.

3 Software architecture

This section provides a comprehensive overview of PaleoSTeHM, detailing its foundational model implementation (Sect. 3.1), the basic architecture for a typical PaleoSTeHM experiment (Sect. 3.2), and the development of PaleoSTeHM modules (Sects. 3.3, 3.4, and 3.5).

3.1 Model implementation

PaleoSTeHM is designed to be a functionally extensible and high-performing toolkit for modeling paleo-data. It is fully open-source and developed under a four-layer structure to maintain a flexible and generic design that is agile to future development (Fig. 1). The four layers, shown from bottom to top in Fig. 1) are (L1) computing platforms, (L2) machine learning platforms, (L3) PaleoSTeHM modules, and (L4) PaleoSTeHM users. At the fundamental level, PaleoSTeHM utilizes computational power from various platforms (L1), such as clouds, clusters, and high-performance computing systems, to ensure scalability and flexibility for diverse applications. Built upon L1, L2 employs Python as the user interface language, leveraging PyTorch (Paszke et al., 2019) and Pyro (Bingham et al., 2019) as its high-performance machine learning platforms. The fast-evolving ecosystem of these popular machine learning platforms enables PaleoSTeHM to support probabilistic programming, auto-differentiation, GPU acceleration, and state-of-the-art optimization algorithms, making it highly efficient and adaptable to a variety of paleo-environmental statistical tasks.

Figure 1 Schematic illustration of the four-layer structure of PaleoSTeHM. L1 specifies various computing platforms (clouds, clusters, and HPC); L2 comprises machine learning platforms (Pyro (Bingham et al., 2019) and PyTorch (Paszke et al., 2019)); L3 includes PaleoSTeHM modules (modeling choices, GP kernel, and analysis choices; see Fig. 2); and L4 consists of the user layer, facilitating interaction with external packages and tools for practical applications and scientific inquiries.

The core toolkit and development reside in L3, which comprises modules that integrate existing machine learning capabilities from L2. This layer includes three primary components: (1) the modeling choices module, which provides options for data-, process-, and parameter-level modeling (Sect. 3.3); (2) the Gaussian process kernel module, a sub-module of the modeling choices module that supports kernel construction using GP priors (Sect. 3.4); and (3) the analysis choices module, which incorporates methods to propagate temporal uncertainty into inference results (i.e., temporal uncertainty treatment) and Bayesian inference (Sect. 3.5). These modules enable flexible and efficient spatiotemporal hierarchical modeling for a wide range of paleo-environmental applications.

We anticipate PaleoSTeHM interacting with external packages and/or tools for practical applications and addressing scientific questions on the PaleoSTeHM user layer (L4, Fig. 1). Here, “External Packages” refers to external Python libraries, which provide various pre-processing and post-processing data functions. For example, in PaleoSTeHM tutorials (see Sect. 4), we use SciPy (Virtanen et al., 2020) for interpolation and Matplotlib (Hunter, 2007) for visualization. “Tools” represents frameworks and services adapted by other developers to integrate PaleoSTeHM capabilities into their toolkits (e.g., Framework for Assessing Changes To Sea-level (FACTS); Kopp et al., 2023). Such plug-in implementations will make it easy for users drawn from any of the PaleoSTeHM categories to use, extend, or contribute to core capabilities for various scientific applications.

3.2 PaleoSTeHM modeling workflow

Constructing and optimizing a hierarchical model within PaleoSTeHM involves a workflow consisting of five sequential selection steps (outlined in Fig. 2), with a focus on modeling and analysis choices in layer L3 as depicted in Fig. 1. Typical PaleoSTeHM experiment steps include (1) selecting data-level models for paleo-environmental data, (2) choosing an appropriate process-level model to describe the latent process, (3) defining prior distributions for each model parameter, (4) selecting a temporal uncertainty treatment method, and (5) choosing a Bayesian inference method (Fig. 2). These five steps reflect core functionalities of PaleoSTeHM modules (layer L3 shown in Fig. 1). To support the effective selection of modeling and analytical choices provided by PaleoSTeHM for various paleo-environmental applications, the fundamental theories and example applications for each modeling choice are introduced in Sect. 3.3.

Figure 2 A schematic illustration of the PaleoSTeHM modeling workflow, providing more detailed information about layer L3 in Fig. 1. The large numbered boxes represent five steps to build a hierarchical model, and it should be noted that the data-level model is specified within each process-level model in PaleoSTeHM v1.0. The smaller boxes indicate different modeling choices within each step. Gray boxes denote available choices that apply to other gray boxes in different steps. Red and purple boxes represent a specific data-level model and temporal uncertainty treatment method corresponding to a specific process-level model (e.g., temporally linear and Gaussian process models), as indicated by colored arrows. The dashed gray box (“Physical model”) highlights that no specific physical model is implemented in PaleoSTeHM. Instead, PaleoSTeHM utilizes outputs from other physical models (see Sect. 3.3.2).

3.3 Modeling choices module

As mentioned above, spatiotemporal hierarchical modeling experiments begin with selecting an appropriate modeling choice for a specific problem. This module offers a variety of commonly used temporal and spatiotemporal modeling options for paleo-environmental studies (Fig. 2). While this paper does not include a dedicated section for parameter-level modeling, the integration of Pyro (Bingham et al., 2019) and PyTorch (Paszke et al., 2019) enables users to define prior probabilities for data- and process-level model parameters using a wide range of commonly used probability distributions. This functionality allows users to customize priors as needed for their specific modeling requirements.

3.3.1 Data-level modeling

The data level of a hierarchical statistical model characterizes the relationship between true (unobserved) target signals and uncertain observations due to multiple error sources. For example, in reconstructing past sea-level changes, the data level addresses uncertainties arising from elevation measurements, indicative range, and leveling errors (Khan et al., 2017). Additionally, proxy data are often subject to inherent temporal uncertainties stemming from various geochronological methods (e.g., radiocarbon dating; Reimer et al., 2020; Heaton et al., 2020). This relationship between observed data and latent process can be formally expressed as

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{y}_i=f(x_i,t_i)+{\mathit{ϵ}}_i^y,\text{(6)}& {\displaystyle }& {\displaystyle }{t}_i={\widehat{t}}_i+{\mathit{ϵ}}_i^t,\end{array}$$

where y_i is the observed data, x_i is the noise-free spatial location of ith observation, t_i is its true age, ${\widehat{t}}_i$ is the mean observational age, and ${\mathit{ϵ}}_i^t$ and ${\mathit{ϵ}}_i^y$ are uncertainties in the age measurement and target signal reconstruction. For paleo-environmental studies, a commonly made assumption is that both ${\mathit{ϵ}}_i^t$ and ${\mathit{ϵ}}_i^y$ are multivariate normally distributed with zero mean and heteroscedastic covariance, so ϵ^y can be expressed as

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{\mathit{ϵ}}^y\sim \mathcal{N}(\mathrm{0},{\mathrm{\Sigma }}_y),\text{(8)}& {\displaystyle }& {\displaystyle }{\mathrm{\Sigma }}_y=\left(\begin{array}{cccc}\text{var}\left(y_{\mathrm{1}}\right)& \text{cov}(y_{\mathrm{1}},y_{\mathrm{2}})& \mathrm{\cdots }& \text{cov}(y_{\mathrm{1}},y_n)\\ \text{cov}(y_{\mathrm{2}},y_{\mathrm{1}})& \text{var}\left(y_{\mathrm{2}}\right)& \mathrm{\cdots }& \text{cov}(y_{\mathrm{2}},y_n)\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \text{cov}(y_n,y_{\mathrm{1}})& \text{cov}(y_n,y_{\mathrm{2}})& \mathrm{\cdots }& \text{var}\left(y_n\right)\end{array}\right),\end{array}$$

where n indicates the number of observations available; var(⋅) represents the variance of specific data; and $\text{cov}(\cdot ,\cdot )$ stands for covariance between two data points, which is often assumed to be 0 when all data are assumed to be independently distributed. However, in practice, strong correlations in paleo-environmental data can emerge from shared processes or dependencies, such as sedimentary records from the same core or data dated using age–depth modeling techniques, where shared depositional history introduces correlated uncertainties (Cahill et al., 2015; Blaauw, 2010). Ignoring these correlations can lead to biased estimates and reduced model reliability. Adapting the likelihood structure to account for covariance, for example, by using a structured covariance model from age–depth modeling, allows more accurate and robust inference (Cahill et al., 2015).

In PaleoSTeHM v1.0, the data-level model is specified within each process-level model, which is assumed to be normally and independently distributed (Fig. 2). For illustrative purposes, PaleoSTeHM v1.0 also includes an implementation of uniform likelihood together with a temporally linear model (see Fig. 2 and Sect. 4.1). For specific problems requiring different likelihood structures, users can replace the likelihood sampling code (a probabilistic random-sampling operation in Pyro) to utilize most of the standard probability distributions supported by Pyro, such as multivariate normal distributions with covariance structures mentioned above.

3.3.2 Process-level modeling

The process level is a hierarchical layer where the variability in the paleo-environment signal is modeled and, in certain cases, decomposed. The process level reflects a scientific understanding of environmental change processes. PaleoSTeHM v1.0 offers multiple process-level models for temporal or spatiotemporal data analysis.

Temporally linear models. Starting with temporal data analysis, probably the most straightforward method for estimating linear trends and the average rate of paleo-environmental change is to fit a linear model to the observed data over time (i.e., straight-line model). For example, Engelhart et al. (2009) and Islam et al. (2021) applied linear regression to discrete paleo-environment data to estimate, respectively, the average rate of sea-level, rainfall, and temperature change during specific time intervals. Over those periods, the observations were assumed to be well represented by a linear trend. A temporally linear model can be expressed as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }f\left(t\right)=\mathit{\alpha }+\mathit{\beta }t,\end{array}$$

where f(t) is the modeled true RSL, β is the constant rate of change in paleo-environmental variable, and α is the intercept (Ashe et al., 2019).

Change-point models. Change-point models describe a single time series by partitioning it into distinct, contiguous segments, each characterized by a linear trend over time (Carlin et al., 1992). These models are widely used to identify the timing of abrupt changes in past climate conditions. For instance, Caesar et al. (2021) and Kemp et al. (2015), respectively, employed change-point models to determine the onset of reduced strength in the Atlantic Meridional Overturning Circulation and the commencement of modern sea-level rise in Connecticut. With m change points, the change-point model can be written as

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }f\left(t\right)=\left\{\begin{array}{ll}{\mathit{\alpha }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{1}}(t-{\mathit{\gamma }}_{\mathrm{1}}),& \text{when }t<{\mathit{\gamma }}_{\mathrm{1}}\\ {\mathit{\alpha }}_{j-\mathrm{1}}+{\mathit{\beta }}_j(t-{\mathit{\gamma }}_{j-\mathrm{1}}),& \text{when }{\mathit{\gamma }}_{j-\mathrm{1}}<t<{\mathit{\gamma }}_j\\ {\mathit{\alpha }}_m+{\mathit{\beta }}_{m-\mathrm{1}}(t-{\mathit{\gamma }}_m),& \text{when }{\mathit{\gamma }}_m\le t,\end{array}\right.\end{array}$$

where γ_k represents a change point, α_k denotes the expected value of RSL at that change point, and β_j indicates the rate of RSL change for each of the m+1 segments. This model incorporates a continuity constraint ensuring that α_k equals α_k−1 plus the product of β_k−1 and the difference between γ_k and γ_k−1. In PaleoSTeHM, the change-point model is implemented to allow users to specify any number of change points (i.e., m in Eq. 10) in the model.

Gaussian process models. GP modeling is a non-parametric Bayesian approach that has been frequently used to infer temporal (or spatiotemporal) variation in paleo-environmental change, including magnitude and rate (Ashe et al., 2019). In models with GP priors, the relationships among any set of points (e.g., over time or across both space and time) are described by a multivariate normal distribution, fully characterized by a mean function and a covariance function (or kernel). Unlike parametric models, such as linear or change-point models used for spatiotemporal analysis, GP models offer greater flexibility because the shape of the curve is determined by the covariance matrix, which reflects the relationship between data points and is inferred directly from the data, rather than being restricted by a predefined functional form.

GP models have gained considerable traction in paleo-environmental science, largely owing to their proficiency in extracting meaningful insights from relatively small datasets. They utilize a non-parametric framework to interpret intricate data patterns effectively. For example, Kay et al. (2021) utilized a GP model to assess herbivore richness for different latitudes in Argentina. Apart from that, Walker et al. (2021) estimated the trend and rate of RSL change across the US Atlantic coast with a GP model. A spatiotemporal GP model, which is defined by its mean function, μ(X), and covariance function (i.e., kernel) $K(X,X^{\prime })$, can be expressed as

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }f\left(X\right)\sim \text{GP}\left(\mathit{\mu }\right(X),K(X,X^{\prime }\left)\right),\end{array}$$

where X indicates spatiotemporal location. A popular choice for many paleo-environmental studies is using the zero-mean function, indicating μ(X)=0 everywhere. In this case, the predictions are only determined by covariance function $K(X,X^{\prime })$, which defines prior expectations about how information is shared between points in different time and space, which typically decays as the time and space differences increase (Rasmussen and Williams, 2006).

Constructing the covariance function is a pivotal and challenging step in a GP model, as it significantly influences the outcome of the inference results. However, justifying the form of the covariance function in Gaussian processes for paleo-environmental studies can be challenging because the processes being modeled are influenced by a wide range of spatial and temporal dependencies, many of which are complex, non-stationary, and not well understood (Tingley et al., 2012; Stein, 2012, 2005a). PaleoSTeHM addresses this by incorporating a “GP kernel” module under the modeling choices module, designed to offer more flexibility and customization extendability. This module provides a user-friendly platform for creating and managing GP kernels, streamlining the process of model construction, and enhancing the adaptability of the analysis to diverse problems. For paleo-environmental applications, multiple choices of building kernels have been adopted in various studies (e.g. Walker et al., 2021; Hay et al., 2015; Kopp et al., 2016, 2014, 2009), and some examples are shown in Sect. 4.2.

Physical models. A physics-based model simulates real-world changes with predictive capabilities anchored in the causal mechanisms delineated by the laws of physics (Saltzman, 2001; Farrell and Clark, 1976). Comparatively, statistical models mostly depend on data-driven correlations, often overlooking fundamental physical principles (e.g., mass or energy conservation). Examples in paleo-environment research include using global circulation models to understand the response of the climate system to different climate forcings (Kageyama et al., 2018) and employing ice sheet dynamic models to quantify past ice sheet response to climate change (DeConto and Pollard, 2016; Tarasov et al., 2012). In the realm of modeling paleo-sea-level change, the GIA model is a widely adopted tool to characterize sea-level changes driven by the gravitational, rotational, and deformational (GRD) effects resulting from the redistribution of ice and water mass (e.g. Lin et al., 2023a; Whitehouse, 2018). The predictive power of such a model is contingent upon underlying formulation and core physical parameters (Peltier et al., 2015; Kendall et al., 2005; Peltier, 2004), such as the history of ice sheet fluctuations and the rheological properties of the Earth's interior for a GIA model (Lin and Yousefi, 2025; Austermann et al., 2013). Validating the physical model against observational data should allow a more accurate representation of spatially dependent patterns of sea-level change, including those linked to sea-level fingerprints (Lin et al., 2021; Mitrovica et al., 2001), in stark contrast to statistical models that might merely presume correlation diminishes with distance (Walker et al., 2021).

Although PaleoSTeHM does not include a specific type of physics-based model (Fig. 2), it offers multiple options to incorporate physical model outputs into final estimates (see examples in Sect. 4.2). Users can use PaleoSTeHM to probabilistically calibrate physical model ensembles conditioned upon observational data. For instance, latent paleo-environmental processes can be modeled as a combination of physical model ensembles conditioned on different physical parameter combinations, using a Dirichlet distribution prior. PaleoSTeHM also supports using a physical model as a mean function in a GP model. In this context, the GP covariance function essentially models the residuals – those processes not captured by the physical model – between observations and the predicted mean function. Additionally, PaleoSTeHM facilitates the construction of sampling covariance functions derived from a physical model ensemble, further enhancing its utility in model integration and assessment (Hay et al., 2015). All of these capabilities are demonstrated in Sect. 4.2, with accompanying source code provided on the PaleoSTeHM GitHub page (see Code and data availability).

3.4 Gaussian process kernel module

The GP kernel module in PaleoSTeHM is a cornerstone for modeling spatial and temporal variations in paleo-environmental data based on GP priors (Figs. 1 and 2). It includes a variety of widely used kernels as described in Rasmussen and Williams (2006). In paleo-environmental studies, examples of kernel applications include the linear (or dot-product) kernel (Khan et al., 2017), radial basis function kernel (Cahill et al., 2015), rational quadratic kernel (Turner et al., 2023; Hay et al., 2015), Matérn kernel (Walker et al., 2021; Kopp et al., 2016), and periodic kernel (Meltzner et al., 2017). These kernels characterize features such as stationarity, isotropy, smoothness, and periodicity in Gaussian processes (Ashe et al., 2019). Detailed kernel information is given in Table 2.

Table 2 Summary of Gaussian process kernels in PaleoSTeHM, based on Stein (2012) and Rasmussen and Williams (2006).

^a All GP kernels can calculate temporal covariance, except the spatial white noise kernel. ^b X represents spatiotemporal location, incorporating both the age and spatial coordinates of the data; t denotes the age of the sample; and x indicates the spatial coordinates. σ² is variance; ℓ is the positive-characteristic length-scale parameter; α is the scale mixture parameter, when α→∞ and the rational quadratic kernel is equivalent to the radial basis function kernel; γ is the offset or shift parameter, adjusting the baseline level of the kernel's output; p is the periodicity parameter for the periodic kernel, defining the cycle length of repeating patterns; and λ represents the degree of the polynomial, an integer determining the complexity of the model for the polynomial kernel. ^c $\text{Cov}\left(m\right(X),m(X^{\prime }\left)\right)$ indicates the sampling covariance between outputs at different spatiotemporal points, derived from deterministic models under varying physical parameter assumptions. Here, m(X) denotes the output at a specific location and time from a suite of physical models assuming different parameters.

Each kernel possesses unique characteristics and necessitates specific parameters (Table 2). For instance, the linear kernel produces linear trends identical to a temporally linear model, suitable for modeling signals with long temporal length scales (e.g., tectonic and GIA in Common Era and future sea-level modeling; Kopp et al., 2016, 2014). The radial basis kernel and the Matérn family of kernels are highly generalizable and allow specification of the degree of differentiability (Table 2), making them suitable for representing physical processes with different levels of smoothness. For example, the GRD effects related to GIA are spatiotemporally smooth, while sediment-compaction-induced sea-level rise can be much more localized and rough (i.e., less differentiable; Kopp et al., 2016; Mitrovica et al., 2011).

In the GP kernel module of PaleoSTeHM v1.0, all kernels are designed for process-level modeling to capture temporal and/or spatial correlations, except for the temporal and spatial white noise kernels, which account for additional measurement errors or unstructured variability by introducing serially uncorrelated uncertainty at the data level (Eq. 5). Apart from the linear kernel, all included kernels are stationary and isotropic (Table 1). To enhance kernel construction flexibility, PaleoSTeHM supports combining different kernels, either additively, multiplicatively, or both. Additive combinations capture independent contributions from distinct processes, such as long-term trends or periodic variations, treating them as separate effects. In contrast, multiplicative combinations create interactions between processes, resulting in more structured patterns. For example, multiplying a periodic kernel with a linear kernel produces a periodic variation with an amplitude that increases or decreases linearly over time, effectively modeling phenomena where seasonal patterns intensify or diminish progressively (Görtler et al., 2019).

Designed for spatiotemporal data analysis, all GP kernels in PaleoSTeHM support temporal data (represented as a one-dimensional (1D) vector), and most of the kernels support spatial data (represented as a two-dimensional matrix including latitude and longitude; see Table 2). Temporal kernel correlations are calculated using the one-dimensional Euclidean distance between time points, while spatial kernel correlations are derived from the one-dimensional geographical radial distance, calculated based on the spherical distance between pairs of longitude and latitude under the assumption of a purely spherical Earth geometry. Users can choose to build a temporal or spatial kernel by switching a parameter in each kernel function.

3.5 Analysis choices module

To accommodate diverse computational resources and varying requirements for the trade-off between modeling robustness and computational demands, the analysis choices module offers multiple methods for Bayesian inference of model parameters as defined in the modeling choices module (Fig. 2). This flexibility ensures users can optimize their analyses based on available technology and specific modeling needs. Unlike deterministic methods (e.g., least-squares), which have been extensively implemented in other studies (e.g. Crichton et al., 2023; Lin et al., 2021), PaleoSTeHM focuses on developing Bayesian probabilistic approaches that more effectively manage the inherent uncertainties associated with paleo-data.

3.5.1 Fully Bayesian analysis

A fully Bayesian analysis requires assigning prior probability distributions to all model parameters, allowing them to take on a range of values, potentially with different probabilities. These priors can either incorporate informative prior knowledge or remain uninformative and vague. Since the posterior distribution is shaped by both the priors and the likelihood of the observed data, it often becomes complex and analytically intractable. Markov chain Monte Carlo (MCMC) methods are crucial in this case, as they enable the efficient exploration and approximation of the posterior distribution. PaleoSTeHM supports two advanced MCMC samplers, Hamiltonian Monte Carlo (HMC; Neal, 2011) and the No-U-Turn sampler (NUTS; Hoffman and Gelman, 2014), which provide more efficient sampling performance than traditional Metropolis–Hastings MCMC (Hastings, 1970).

HMC significantly improves sampling efficiency over traditional Metropolis–Hastings MCMC by leveraging gradients of the probability distribution to guide the sampling process, which involves generating random samples from the underlying latent probability distribution. This method reduces autocorrelation (the correlation between successive samples in a Markov chain, indicating how dependent the current sample is on previous ones), thereby increasing the effective sample size (the number of independent samples, accounting for autocorrelation) per iteration (a single step in the sampling process where the algorithm generates a new sample) and enabling faster convergence. Building on HMC, NUTS further enhances efficiency by automatically adapting the path length (the distance traversed in parameter space during a single Hamiltonian trajectory) and managing the step size (the distance traveled in parameter space at each leapfrog step during Hamiltonian dynamics; Bingham et al., 2019). NUTS eliminates the need for manual tuning of these parameters, facilitating more effective exploration of complex, high-dimensional posterior distributions commonly encountered in Bayesian analysis.

Compared to other analysis choices such as empirical Bayesian models or variational Bayesian models (Table 1), a fully Bayesian model offers a more comprehensive estimation of the uncertainties associated with model parameters (Piecuch et al., 2017). It also offers a direct framework for sample age measurement uncertainty in an errors-in-variables (EIV) manner (Dey et al., 2000). However, the nature of MCMC-based samplers means they are computationally more demanding. Particularly within the EIV framework, where the number of sampling parameters increases linearly with data size, this leads to a polynomial increase in the computational power required (Belloni and Chernozhukov, 2009), which can be significant and unaffordable when dealing with large datasets or complex models.

3.5.2 Empirical Bayesian analysis

Unlike fully Bayesian analysis, which requires full probability distributions for prior and posterior, empirical Bayesian analysis offers a practical alternative. This approach approximates a fully Bayesian treatment where parameters at the highest level of the hierarchy are fixed at their most likely values rather than being integrated out. This optimization is typically achieved using the maximum likelihood estimate, leading to a posterior distribution that is conditional on the data and these optimized parameters:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }p\left(f\right|y,{\widehat{\mathit{\theta }}}_s,{\widehat{\mathit{\theta }}}_d)\propto p(y|f,{\widehat{\mathit{\theta }}}_d)p\left(f\right|{\widehat{\mathit{\theta }}}_s).\end{array}$$

Here, the posterior probability of the latent processes f is inferred, assuming that the hyperparameters at the data and process levels (${\widehat{\mathit{\theta }}}_d$ and ${\widehat{\mathit{\theta }}}_s$) are known and fixed. While the existing code base allows explicit bounds to be set on hyperparameters for the maximum likelihood estimate (e.g. Ashe et al., 2019; Kopp et al., 2016), it does not provide for an explicit prior distribution for the parameters. In other words, it only supports uniformly distributed prior information, limiting the ability to incorporate informative prior knowledge. By leveraging Pyro's variational inference capabilities (details in Sect. 3.5.3), PaleoSTeHM enables users not only to optimize hyperparameters using their maximum likelihood estimate but also to define many commonly used distributions for each prior model parameter explicitly. This allows optimization to be conducted in a maximum a posteriori probability estimation manner, assuming the variational distribution is a Dirac delta function. In PaleoSTeHM, by default, the optimization is achieved using Adam, a stochastic optimizer (Kingma and Ba, 2014). While empirical Bayesian analysis generally requires fewer computational resources than fully Bayesian methods, it is important to note that, assuming hyperparameters at the data and process levels are known and fixed, it may lead to substantial underestimation in the inference uncertainty (Piecuch et al., 2017).

3.5.3 Variational Bayesian analysis

Considering the computational expense required to perform MCMC in fully Bayesian analysis and the limitations of empirical Bayesian methods that fail to account for the uncertainty of hyperparameters, PaleoSTeHM also supports variational Bayesian analysis, which emerges as an efficient intermediary. Rather than directly sampling from the posterior distribution through MCMC, variational Bayesian methods aim to approximate the true posterior probability distribution ($p(f,{\mathit{\theta }}_s,{\mathit{\theta }}_d|y)$) with a simpler, parametric probability distribution (q(f|ϕ)). Thus, Bayesian inference is transformed from a sampling challenge into an optimization problem – known as variational inference – requiring significantly fewer computational resources while facilitating uncertainty estimation (Wingate and Weber, 2013).

In PaleoSTeHM, variational Bayesian analysis is achieved by optimizing the variational parameters ϕ to minimize the Kullback–Leibler (KL) divergence, a metric to effectively measure the difference between two distributions:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\mathit{\varphi }=\arg \underset{\mathit{\varphi }}{min}\text{KL}\left[q(f,{\mathit{\theta }}_s,{\mathit{\theta }}_d|\mathit{\varphi }\left)\right|\left|p\right(f,{\mathit{\theta }}_s,{\mathit{\theta }}_d\left|y\right)\right].\end{array}$$

For more details above KL divergence, readers can refer to Blei et al. (2017). Adam facilitates this minimization, and the variational distribution for PaleoSTeHM is a normal distribution by default. In contrast to MCMC-based fully Bayesian analysis, which often requires computational power that increases polynomially with the number of data points, the optimization-driven approach of variational Bayesian analysis generally scales linearly (Ko et al., 2024; Hoffman and Blei, 2015). Consequently, variational methods can handle larger datasets more effectively, making them suitable for large-scale problems prohibitively for full Bayesian analysis.

3.5.4 Incorporation of temporal uncertainty

PaleoSTeHM provides two methods to incorporate temporal uncertainty into final estimations. The first method uses EIV framework (Cahill et al., 2015; Dey et al., 2000), which directly incorporates temporal uncertainty through MCMC sampling of the distribution. The second approach adopts the noisy-input framework (McHutchon and Rasmussen, 2011), which applies a first-order Taylor series approximation – a linear expansion around each input point – to account for errors in the independent variable, time, thereby translating these into equivalent errors in the dependent variable:

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }f(x_i,t_i)\approx f(x_i,{\widehat{t}}_i)+{\mathit{ϵ}}_i^t{\displaystyle \frac{\partial f(x_i,{\widehat{t}}_i)}{\partial t}}.\end{array}$$

Here, ${\widehat{t}}_i$ and ${\mathit{ϵ}}_i^t$ are the same as in Eq. (6), standing for mean observational age and age uncertainty, respectively. The integration of temporal uncertainty within PaleoSTeHM is executed alongside each process-level model (Fig. 2). All process-level models are implemented using an EIV framework, while, for the GP models, both EIV and noisy-input frameworks are available (Fig. 2).

3.6 Model validation

After implementing and optimizing a hierarchical model in PaleoSTeHM, it is essential to perform a model validation step to further ensure the robustness and reliability of the trained model. This process involves evaluating how well the model fits the observed data, assessing its predictive accuracy, and diagnosing potential issues such as convergence problems. PaleoSTeHM includes a range of techniques for model validation, such as residual analysis, posterior predictive checks, MCMC convergence diagnostics (e.g., effective sample size and Gelman–Rubin statistic; Gelman and Rubin, 1992), visual inspections (e.g., optimization trace plots, true vs. predicted plot), and simulation validation and cross-validation methods. These tools allow users to critically examine the model's assumptions, quantify uncertainties, and compare competing models to select the most appropriate one for their specific paleo-environmental application.

To complement these validation techniques, we demonstrate their application in various case studies presented in Sect. 4. Each case study incorporates specific model validation methods tailored to the modeling and analysis choices used. For example, prior and posterior predictive checks are employed to evaluate the performance of optimized models (Sect. 4.1.1); residual plots, weighted mean squared error (wMSE), and cross-validation are used to assess the performance of different process-level models (Sects. 4.1.2 and 4.2); and effective sample size and the Gelman–Rubin statistic ensure good model convergence for MCMC-based analyses (Sect. 4.1.3). Detailed implementations and usage of these validation methods for various PaleoSTeHM experiments, including those mentioned above but not covered in detail in the following sections, are available in the PaleoSTeHM tutorials (see Code and data availability).

4 Case studies

This section presents illustrative case studies using a tutorial format to demonstrate PaleoSTeHM's usability. All codes and data are accessible and actively managed on the PaleoSTeHM GitHub page (https://github.com/radical-collaboration/PaleoSTeHM, last access: 1 April 2025). The case studies include the following:

1. reconstruction of temporal sea-level changes using coral reef data from the Great Barrier Reef with different data-level models (Sect. 4.1.1),

2. reconstruction of temporal sea-level changes using salt marsh data from New Jersey with different process-level models (Sect. 4.1.2),

3. reconstruction of temporal sea-level changes using salt marsh data from North Carolina with different Bayesian inference methods (Sect. 4.1.3),

4. reconstruction of spatiotemporal sea-level changes using various geological proxies from the US Atlantic coast with different process-level models (Sect. 4.2).

Although these examples focus on modeling paleo-sea level, additional tutorials are available for analyzing other paleo-environmental data, such as ocean temperature anomalies and concentration of carbon dioxide.

The prior and posterior distributions and analysis choice for each model are provided in Table A1. It should be noted this section only briefly describes the modeling results; for a more systematic analysis of paleo-environmental modeling results based on different statistical techniques, the user can refer to Ashe et al. (2019), PAGES2k Consortium (2019), and Tingley et al. (2012).

4.1 Time series analysis

4.1.1 Data-level modeling

In this section, we examine the impact of the data-level model on inference results. Although numerous paleo-environmental applications commonly assume that proxy reconstruction uncertainties are normally distributed (Ashe et al., 2019; Khan et al., 2019; Tingley et al., 2012), certain types of proxies may exhibit different forms of uncertainty. For instance, coral reef sea-level proxies indicate past sea-level changes through a quantifiable relationship between the coral's living-habitat depth and the concurrent sea level (Hibbert et al., 2016). Representations of coral living-habitat depth uncertainties are often modeled using either a normal distribution (e.g. Khan et al., 2019) or a uniform distribution (Lin et al., 2021). To illustrate such data-level impact on inference results, we apply a temporally linear model within an EIV framework to coral reef data from the Great Barrier Reef (Yokoyama et al., 2018) using two alternative data-level models. The first can be expressed as

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{\mathit{ϵ}}_{\mathrm{1}}^y\sim U({\mathit{\tau }}_{\mathrm{l}}-{\mathit{\omega }}_{\mathrm{1}},{\mathit{\tau }}_u+{\mathit{\omega }}_{\mathrm{1}}),\end{array}$$

where U indicates a uniform distribution between lower and upper ranges defined by specific coral species (τ_l and τ_u) and an additional noise, defined by hyperparameter ω₁, which follows a prior distribution of

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{\mathit{\omega }}_{\mathrm{1}}\sim U({\mathrm{10}}^{-\mathrm{4}},\mathrm{10}).\end{array}$$

The second data-level model can be represented as

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{\mathit{ϵ}}_{\mathrm{2}}^y\sim N({\mathit{\mu }}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}),\end{array}$$

where N indicates a normal distribution with mean μ₂ and a standard deviation σ₂, both of which are determined by specific coral species, and ω₂ is an additional noise hyperparameter, following the same prior distribution as ω₁. The same prior distributions for each parameter are used for both data-level models, which are represented as non-informative uniform distributions. The characteristics of these non-informative priors are evident in the prior predictive check (Fig. A1), which reveals a wide and flat spread of predictions, reflecting the absence of observational influence at this stage.

For both models, the posterior distribution is determined by 11 000 posterior samples drawn from a NUTS sampler, with the first 1000 samples discarded as burn-in steps. The posterior predictive checks (Fig. A1) illustrate that the posterior predictions for both models align closely with the observed data, suggesting successful model convergence. It can be seen in Fig. 3 that, although the inference results from different data-level models are overall similar, there are still some noticeable differences in the inferred sea-level change trend and rate. Uniform and normal likelihoods yield average sea-level rates of 5.91 mm yr⁻¹ (4.45–7.38 mm yr⁻¹; 90 % credible interval; CI) and 6.29 mm yr⁻¹ (4.81–7.73 mm yr⁻¹), respectively. These likelihood assumptions also produce considerably different additional noise parameter distributions. Therefore, users should select an appropriate data-level model to better represent the specific characteristics of different paleo-environmental data.

Figure 3 The impact of data-level models on temporal sea-level change inference at the Great Barrier Reef. Posterior probability density functions of sea-level change rate (a), intercept (b), and standard deviation of additional noise (c), assuming either a uniform likelihood (blue) or a normal likelihood (orange). Inferred mean sea-level trends (d) and rates (e) along with a 90 % credible interval, where sea-level data are represented by red boxes, with horizontal range indicating ±2σ age uncertainty and vertical range indicating the reconstructed maximum and minimum sea-level range determined by coral species (i.e., τ_l and τ_u in Eq. 15). A negative RSL value indicates that the local RSL in the Great Barrier Reef was lower than present-day levels, reflecting the significant amount of water stored in continental ice sheets. CI: credible interval. BP: before present.

4.1.2 Process-level modeling

To demonstrate the impact of different process-level models on inferring paleo-sea-level time series, we use the same data-level model together with process-level models introduced in Sect. 3.3.2, employing non-informative priors (Table A1). The sea-level data used here are a near-continuous core record from single cores of salt marsh sediment from Leeds Point (New Jersey) covering the Common Era (Kemp et al., 2013). For this database, a normal-likelihood data-level model is adopted with sea-level reconstruction uncertainties provided by the original study. Here, we test three process-level models: (a) a temporally linear model, (b) a change-point model (assuming three change points; following Ashe et al., 2019), and (c) a Gaussian process model with a radial basis function (RBF) kernel. Posterior distributions for models (a) and (b) were obtained using a variational Bayesian approach, while model (c) employed an empirical Bayesian method.

Figure 4 shows estimated RSL trends and rates of RSL change for each process model. The resulting trends and fit to the data (quantified using wMSE) differ significantly due to the fundamentally different model formulations. The temporally linear model can only estimate an averaged trend and rate of sea-level change and will never predict an accelerated RSL change. Consequently, it exhibits the highest wMSE (2.53) with systematic errors with strong temporal correlations displayed in residual plots (Fig. A2), both of which reflect a poor fit to the observations.

Figure 4 The impact of process-level models on temporal sea-level change inference at New Jersey. Common Era sea-level comparison of a linear model (a, e), a change-point model (assuming three change points; b, f), and a Gaussian process model with an RBF kernel (c, g), where input data are continuous cores. Output includes estimates of RSL (a–d) and rates of RSL change (e–h), which are each shown with mean and 90 % credible intervals. Paleo-sea-level data are modeled here using a normal likelihood. The horizontal and vertical ranges of red boxes indicate ±2σ age and relative sea-level reconstruction uncertainties, respectively. Sea-level data were reconstructed here using a near-continuous record from single cores of salt marsh sediment from Leeds Point (New Jersey; Kemp et al., 2013). wMSE: weighted mean squared error. CE: Common Era.

Comparatively, the change-point model is able to capture a noticeable change in RSL rate from 1.49 mm yr⁻¹ (1.26–1.70 mm yr⁻¹) between −500 CE and 1839 CE (1824–1852 CE) to 3.91 mm yr⁻¹ (3.72–4.10 mm yr⁻¹) after 1839 CE. This added flexibility substantially improves the model's fit to the data, achieving a wMSE of 0.38 and producing a less structured error distribution (Fig. A2). Such flexibility makes the change-point model particularly suitable for identifying the time of emergence in various environmental change contexts (e.g. Walker et al., 2022; Caesar et al., 2021; Lyu et al., 2014).

As a non-parametric approach, the GP model produces continuous distributions of RSL change rates over time, allowing the estimation of multiple inflection points (Walker et al., 2022). This flexibility results in the lowest wMSE (0.31), alongside minimal temporal structure in the residuals (Fig. A2), indicating the best overall fit to the observations. However, the infinite differentiability of the RBF kernel can lead to overly smooth predictions in time series analysis, potentially oversmoothing sharp changes that are critical in many environmental contexts, such as abrupt sea-level rise (Lin et al., 2021), ocean circulation slowdowns (Caesar et al., 2018), and extreme events like heavy rainfall (Stein, 2012). Alternative kernels (e.g., Matérn kernels) can provide alternative levels of differentiability.

4.1.3 Analysis choices

Using similar near-continuous core data from Sand Point, North Carolina (Kemp et al., 2011), we illustrate the effects of analysis choices on RSL inference. Here, we only use a subset of the original data to better demonstrate the difference between various analysis choices. The adopted data- and process-level models employ a normal likelihood with a GP model using an RBF kernel (Table A1). The hyperparameters are sampled using empirical, fully Bayesian. and variational Bayesian methods. For the fully Bayesian method, the posterior distribution is determined by 5500 posterior samples drawn from a NUTS sampler, with the first 200 samples discarded as burn-in steps. For the empirical and variational Bayesian methods, the hyperparameters were optimized using the Adam optimizer over 1000 iterations (Kingma and Ba, 2014). The run times of each implementation are reported on a 2023 MacBook Pro with an Apple M2 Pro chip.

For MCMC-based fully Bayesian analysis, PaleoSTeHM employs the Gelman–Rubin statistic (Gelman and Rubin, 1992) to verify that the Markov chains have converged to a stationary phase, indicating good convergence. Additionally, the effective sample size is used to assess the amount of information retained, accounting for the correlation in the sequence (Bürkner, 2017). Typically, a Gelman–Rubin statistic of less than 1.1 and an effective sample size greater than 1000 suggest reliable sampling of the posterior distribution. In this case, the analysis meets these criteria with a Gelman–Rubin statistic of 1.0 and an effective sample size exceeding 3000. For empirical and variational Bayesian methods, validation is typically conducted through the inspection of optimization trace plots (plots showing how optimization target function improve with each iteration), where successful optimization is characterized by a steadily decreasing loss function and parameter convergence over iterations. These conditions are also satisfied in this analysis, as illustrated in the corresponding optimization trace plots (Figs. A3 and A4).

Figure 5 compares posterior distributions of RSL trend and rate of change and the computational time for each analysis choice. The empirical Bayesian method requires the least computational power, only providing a point estimate of hyperparameters without accounting for their underlying uncertainty. Although more computationally demanding, the fully Bayesian method captures the hyperparameter uncertainties effectively. As an intermediary, the variational Bayesian method requires slightly more computational time compared to the empirical method but can derive a variational posterior distribution that is largely similar to that obtained by the fully Bayesian method through MCMC sampling. In contrast, the point estimate obtained through the empirical Bayesian method falls at the third percentile of the posterior hyperparameter distributions derived from the fully Bayesian method, highlighting a significant bias introduced by the overly simplistic approach.

Figure 5 The impact of analysis choices on temporal sea-level change inference at North Carolina. (a, b) GP model hyperparameter optimization results along with the required computational time in seconds (based on a 2023 MacBook Pro with an Apple M2 Pro chip). (c, d) Common Era sea-level comparison between three analysis choices; the results indicate 90 % credible interval of RSL change trend (c) and rate (d). Paleo-sea-level data are modeled here using a normal likelihood. The horizontal and vertical ranges of red boxes indicate ±2σ age and relative sea-level reconstruction uncertainties, respectively. Sea-level data were reconstructed here using a near-continuous record from single cores of salt marsh sediment from Sand Point (North Carolina, Kemp et al., 2011).

Because of the near-continuous sea-level data with smoothly rising sea-level trend in North Carolina, the inference results from these three methods are similar. However, given that geological sea-level data are often sparsely distributed across both spatial and temporal domains and may be subject to an abrupt change in rate, neglecting the underlying uncertainty of hyperparameters by the empirical Bayesian method can result in a significant underestimation of the final inference uncertainty compared with the fully Bayesian method.

4.2 Spatiotemporal analysis

Spatiotemporal analysis presents a common challenge in paleo-environmental studies, such as reconstructing continuous spatiotemporal signals from sparse and noisy data. To address this, PaleoSTeHM provides a range of approaches, spanning purely statistical to purely physical methods. Here, we present an illustrative example to recover the spatiotemporal RSL pattern and its associated uncertainty. This analysis utilizes a sea-level database containing 1043 proxy records spanning from 11 ka to the present, compiled by Ashe et al. (2019) from previous studies (Kemp et al., 2013, 2014, 2015, 2017a, b; Khan et al., 2017; Engelhart and Horton, 2012). The database includes sea-level proxies such as salt marsh, mangrove, beach rock, and coral. All records were used to train the model except for 51 sea-level data points from New York (Engelhart and Horton, 2012), which were reserved for cross-validation, a technique used to evaluate model performance on unseen data (shown as a gray dot in Fig. 6m–p).

Figure 6 The impact of process-level modeling choices for spatiotemporal sea-level change inference along the US Atlantic coast. Each column represents a different process-level model. (a–d) Mean and 90 % credible or confidence interval for RSL predictions at Florida, also indicated by a red dot in panels (m–p). The horizontal and vertical ranges of red boxes indicate ±2σ age and relative sea-level reconstruction uncertainties, with data from Khan et al. (2017). (b) For model (ii), the GP mean function, determined by the ICE_7G model (Roy and Peltier, 2018; Peltier et al., 2015), is depicted with a black line. (e–h) Normalized residuals, which represent the difference between observed and predicted values normalized by observational uncertainty, for Florida RSL predictions generated by each process-level modeling choice. Dashed red lines represent 0 error. The weighted mean squared error for each model is given above each panel. (i–l) Mean RSL prediction for the year −5500 CE. (m–p) Standard deviation of the RSL prediction for the year −5500 CE, where white dots indicate the locations at which sea-level data were collected and gray dots represent the location of New York (where 51 data points were held for testing purposes).

We demonstrate four process-level models that were used in previous studies: (i) a GP model with a zero mean function and multiple isotropic kernels (Ashe et al., 2019), (ii) a GP model with the mean function determined by a GIA model and multiple isotropic kernels (Walker et al., 2021; Kopp et al., 2016), (iii) a GP model with a zero mean function and a sampling covariance kernel determined by a GIA model ensemble (Kopp et al., 2009), and (iv) a purely GIA model ensemble (Lin et al., 2023a). All models assume a data-level model with a normal likelihood determined by RSL reconstruction uncertainty and an additional white noise term. For this analysis, we implement a noisy-input framework to address temporal uncertainty and use the empirical Bayesian method to optimize hyperparameters for models (i), (ii), and (iii), while model (iv) is optimized through the variational Bayesian method (Table A1).

For model (i), we follow the kernel structure as in Ashe et al. (2019), which can be expressed as

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }f\left(X\right)\sim \text{GP}(\mathrm{0},K_{\mathrm{1}}(X,X^{\prime }\left)\right),\text{(19)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{1}}(X,X^{\prime })=g\left(t\right)+r(x,t)+l(x,t),\end{array}$$

where g(t) represents a spatially uniform covariance function, while r(x,t) and l(x,t) are regional and local varying isotropic covariance functions, respectively. These are characterized by a $\mathrm{3}/\mathrm{2}$ Matérn temporal kernel (Table 2) for g(t) and a product of a $\mathrm{3}/\mathrm{2}$ Matérn temporal kernel and a $\mathrm{1}/\mathrm{2}$ Matérn spatial kernel (Table 2) for r(x,t) and l(x,t), which are distinguished by their prior distributions of hyperparameters.

Similarly, model (ii) can be written as

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }f\left(X\right)\sim \text{GP}\left(\text{GIA}\right(X),K_{\mathrm{2}}(X,X^{\prime }\left)\right),\text{(21)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{2}}(X,X^{\prime })=r(x,t)+l(x,t).\end{array}$$

Here, the mean RSL expectation is determined by RSL prediction from the ICE_7G ice model with the VM5a Earth model (Roy and Peltier, 2018; Peltier et al., 2015), and r(x,t) and l(x,t) are the same as in Eq. (19). We do not include the g(t) kernel here, as the RSL prediction derived from the ICE_7G model – embedded within the GP mean function and rigorously calibrated against comprehensive RSL and geodetic datasets across North America (Roy and Peltier, 2018; Peltier et al., 2015) – is assumed to adequately capture all spatially uniform signals.

Model (iii) can be denoted as

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }f\left(X\right)\sim \text{GP}(\mathrm{0},K_{\mathrm{3}}(X,X^{\prime }\left)\right),\text{(23)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{3}}(X,X^{\prime })=\text{Cov}\left(m\right(X),m(X^{\prime }\left)\right)\cdot \exp (-|t-t^{\prime }|/{\mathit{\tau }}^{\mathrm{2}}).\end{array}$$

Here, Cov(⋅) indicates a sampling covariance function through a physical model ensemble m. In this context, the covariance between data points is not directly determined by their spatiotemporal proximity but instead depends on the variance within the physical model ensemble, conditioned on various combinations of physical parameters. For this (iv), m includes an ensemble of forward GIA models: (a) the ICE_7G ice model with the VM5a Earth model, (b) the PaleoMIST ice model (Gowan et al., 2021) with 71 km lithosphere and 0.3 and 70×10²¹ Pa s upper- and lower-mantle viscosity, and (c) the ANU ice model (Lambeck et al., 2014) with 71 km lithosphere and 1 and 10×10²¹ Pa s upper- and lower-mantle viscosity. To expand the variability in physical model predictions, we create six synthetic GIA model outputs by enlarging or shrinking these three GIA model outputs by 1.5. Therefore, this physical model ensemble consists of predictions from nine models. More details about the physics-based GIA model used here can be found in Lin et al. (2023b). To stabilize the estimate and reduce variability related to finite sample size, we applied a temporal Gaussian taper function to this kernel, controlled by a parameter τ. Following Hay et al. (2015) and Kopp et al. (2009), we set τ to 3000 years.

Lastly, model (iv) can be written as a weighted mean of different physical models:

$$\begin{array}{}\text{(24)}& {\displaystyle }& {\displaystyle }f\left(X\right)=\sum _{n=\mathrm{1}}^N{\mathit{\nu }}_i{\text{GIA}}_i\left(X\right),\text{(25)}& {\displaystyle }& {\displaystyle }\mathit{\nu }\sim \text{Dirichlet}\left({\mathit{\alpha }}_d\right).\end{array}$$

In this model, ν represents the relative weights associated with each GIA model. These probabilities follow a Dirichlet distribution (or multivariate beta distribution) characterized by a concentration parameter α_d. A value greater than 1 for α_d indicates a preference for a more evenly distributed probability across all models. In contrast, a value less than 1 indicates a preference for more concentrated probabilities on fewer models (Lin et al., 2023b). For this experiment, we set α_d according to each GIA model prediction fit to RSL observation (using weighted root mean square as a metric; see Table A1).

A comparison of RSL inference results between different spatiotemporal process-level models is provided in Fig. 6. At the purely statistical end of the process model spectrum, model (i) correlates RSL from various locations and times based solely on their spatiotemporal proximity, a property derived from the adopted isotropic kernels. According to model (i), the RSL change along the US Atlantic coast during the Holocene was dominated by a spatially uniform signal (produced by the spatially uniform kernel, g(t), in Eq. 19; Table A1), which contributed to more than 25 m of RSL rise. In contrast, r(x,t) and l(x,t) only produce up to 5 m of spatially variable RSL signal, resulting in virtually no spatial pattern in the mean RSL prediction of this model. In the temporal domain, multiple studies have demonstrated that GP models like model (i) can accurately recover multi-millennial sea-level variation trends at locations with abundant sea-level observations, such as Florida, as shown in Fig. 6a (Tang et al., 2023; Ashe et al., 2019; Cahill et al., 2015). This is further supported by the residual plot (Fig. 6e), which exhibits a low wMSE (1.1) and minimal temporal structure in the residuals, indicating a good agreement with the observed data.

However, spatial inferences based on isotropic and stationary kernels of model (i) are often considered overly simplistic (Stein, 2005a), partly due to the sparse nature of geological data and the complexity of environmental change mechanisms. As see in Fig. 6m, geological sea-level data are mostly collected across paleo-coastal areas. Therefore, RSL inferences from model (i) are only representative of coastal areas (as opposed to terrestrial or marine areas) and cannot adequately reflect the physical knowledge of paleo-sea-level change (e.g., the RSL uncertainty caused by the existence of the Laurentide Ice Sheet).

Model (ii) uses a deterministic GIA model (ICE_7G ice model with VM5a Earth model; Roy and Peltier, 2018; Peltier et al., 2015) as the GP mean function. By harnessing a physics-based model, model (ii) captures intricate spatial sea-level variation patterns due to the GIA-induced GRD effects (Fig. 4). In this setup, the covariance functions describe residuals between the GIA model and RSL observations (mostly captured by r(x,t) in Eq. 21). Similarly to model (i), RSL predictions for Florida by model (ii) closely align with observations, demonstrating a low wMSE (1.0) and an unstructured residual distribution. In the spatial domain, at −5500 CE, model (ii) suggests the GIA model underestimates ∼10 m RSL at New Jersey (Fig. A5), which may reflect oversimplified physics (e.g., neglecting 3D solid Earth rheology; Austermann et al., 2013), biased sampling of physical parameters (such as poorly constrained ice history), or missing physical processes in the GIA model (e.g., sediment isostatic adjustment; Lin et al., 2023a). Because model (ii) assumes no uncertainty in GIA modeling, the uncertainty quantification here also relies solely on the radial distance from RSL data points (Fig. 6n).

Model (iii) utilizes a kernel constructed by sampling covariances between various forward GIA models based on alternative ice and Earth models. By incorporating relevant physical processes into GP kernel construction, model (iii) effectively captures anisotropic behaviors, non-stationarities, heterogeneities, and teleconnections that are intrinsic to the physical dynamics of RSL change but are challenging to represent with standard classes of covariance functions (Table 2). For instance, the size of the Laurentide Ice Sheet exhibits a positive correlation with RSL around the northern Great Lakes while displaying a negative correlation with RSL in peripheral bulge regions such as New Jersey (Fig. 6k).

Model (iii) also has certain limitations, including the computational cost of thoroughly sampling physical model parameters and the presence of structural errors within the physical models, such as oversimplifications or omitted processes. As shown in the residual plot of Fig. 6g, there is a significant mismatch between RSL observations and model (iii) predictions, reflected in a high wMSE value (17.99) and pronounced temporal structure in the residuals. This poor fit may stem from biased sampling of ice and Earth models and from the reliance on an oversimplified 1D rheology. Furthermore, the model prioritizes fitting regions with denser data distribution, such as the mid-northern US Atlantic coast. For instance, model (iii) provides a good fit to unseen RSL observations in New York (Fig. 7c and g), but this prioritization comes at the expense of accuracy in regions with sparser data, such as Florida, where substantial misfits are observed. Additionally, the posterior mean and standard deviation generated by this method are less directly interpretable compared to those produced by model (iv).

Model (iv) represents the purely physical end of the process-level spectrum and is formulated as a weighted linear combination of physical models, with weights determined by data–model misfits (e.g., wMSE and chi-square misfit; Lin et al., 2021; Li et al., 2020; Lambeck et al., 2014). The mean and uncertainty estimated by this method reflect the parametric uncertainty inherent in a given physical model, allowing direct interpretation of physical parameters, such as deriving posterior distributions of global ice history (Creel et al., 2024). However, this approach is also susceptible to structural errors within the model, similar to those observed in model (iii). The limited sample size of physical parameters – only nine models were used in this analysis – and the model's tendency to prioritize fitting denser sea-level data in mid-northern locations result in uncertainty estimates that appear underestimated and biased, as illustrated by substantial misfits to observations (Fig. 6d and h). Furthermore, the difficulty in directly quantifying certain physical parameters often leads to oversimplified model predictions. For example, the scarcity of direct constraints on ice history (Dalton et al., 2020) reduces the ability of forward GIA models to resolve centennial-scale sea-level variations. This limitation makes these models less effective in capturing centennial-scale variability compared to models (i)–(iii) (Fig. 6).

Due to the dense distribution of sea-level data along the mid-northern US Atlantic coast, all models effectively capture the general trend of RSL variation observed in the withheld data from New York (Fig. 7), as evidenced by low wMSE values and minimal temporal structure in the residual plots. While model (ii) achieves the lowest wMSE and model (iv) exhibits slight bias in its residuals, all models demonstrate comparable overall performance in reconstructing sea-level changes. In contrast, significant variation in model performance is observed for RSL predictions at Florida (Fig. 6a–h), where models (iii) and (iv) show substantial misfit to RSL observations. This misfit stems from the requirement to preserve the overall consistency of the physical systems constrained by the ensemble of physical models.

Figure 7 Performance of different process-level modeling choices on unseen sea-level data from New York. (a–d) Mean and 90 % credible or confidence intervals for RSL predictions at New York (locations shown in Fig. 6m–p as gray dots). The gray boxes represent ±2σ uncertainties in both age and relative sea-level reconstructions, with data from Engelhart et al. (2009). For model (ii), the GP mean function, derived from the ICE_7G model (Roy and Peltier, 2018; Peltier et al., 2015), is shown as a black line (b). (e–h) Normalized residuals, which represent the difference between observed and predicted values normalized by observational uncertainty, for New York RSL predictions generated by each process-level modeling choice. Dashed red lines represent 0 error. The weighted mean squared error for each model is given above the panels.

It is important to note that high-quality and standardized datasets, such as those available for the mid-northern US Atlantic coast, are rare in many paleo-environmental fields, such as deep-sea isotopes or ice core records (Shackleton et al., 2021; Lemieux-Dudon et al., 2010). Consequently, users must carefully evaluate factors such as data availability, computational resources, the need for interpretability, and the level of understanding of underlying physical processes when selecting a process-level model. Generally, physics-based models offer superior interpretability and better extrapolation capabilities to spatiotemporal locations with minimal data, as they are rooted in well-established theoretical frameworks. However, discovering and validating new physical laws can be time-intensive and often computationally demanding. In contrast, machine learning or statistical approaches provide flexibility and computational efficiency but often face challenges in extrapolating nonlinear functions (Xu et al., 2020; Goodfellow et al., 2016), while they require large volumes of training data and rigorous validation to ensure consistency with physical principles. For a more detailed discussion on the integration of physics-based and machine learning models, readers are referred to Lai et al. (2025).

5 Discussion

5.1 Generalization for paleo-environmental problems

This paper focuses on demonstrating the functionality of PaleoSTeHM in paleo-sea-level applications. However, the flexibility of hierarchical models, where any statistical model can be interpreted hierarchically, allows the PaleoSTeHM framework to be readily applicable to a wide range of paleo-environmental problems. By using hierarchical models, transparency is enhanced by distinguishing between modeling assumptions and analytical methods and by separating process variability from observational noise.

The common characteristics of paleo-sea-level datasets, such as sparsity and discreteness, are shared by many other paleo-environmental datasets, including paleo-temperature (PAGES2k Consortium, 2017), past ice sheet thickness (Small et al., 2019), and sediment deposition depth (Wang et al., 2018). As a result, the data- and process-level models introduced in this paper can be readily generalized to these paleo-environmental fields. For example, Tingley et al. (2012) and Stein (2005b) proposed that using a GP model is a reasonable approach to describe latent space–time climate processes, such as annual mean surface temperature anomalies and daily wind speed; Lin et al. (2023a) applied a spatial GP model to recover the spatial pattern of Holocene coral reef depth based on a Holocene coral reef deposition depth database across the Great Barrier Reef (Hinestrosa et al., 2022); and Caesar et al., 2018 implemented a change-point model on multiple proxy datasets to detect significant reductions in the strength of the Atlantic Meridional Overturning Circulation.

Beyond the process-level models featured in PaleoSTeHM v1.0, various approaches have been employed for paleo-environmental analyses. Common techniques for addressing problems in this field include principal component analysis (equivalent to the empirical orthogonal function method when temporal aspects are considered), autoregressive models, and generalized additive models. For instance, Shakun and Carlson (2010) used an empirical orthogonal function approach to detect modes of deglacial temperature variability; Piecuch et al. (2017) adopted a degree-1 autoregressive model to reconstruct sea-level evolution using tide gauge data; and Simpson (2018) and Upton et al. (2023) developed a series of generalized additive models to model paleo-ecology and paleo-sea level, respectively. While the reimplementation of these models in PaleoSTeHM is beyond the scope of this paper, doing so would benefit from the framework's multiple analysis options and its capacity for smooth integration with flexible data- and parameter-level models.

5.2 Future developments

From a scientific perspective, numerous promising directions exist for further development of PaleoSTeHM.

Existing data-level models only support a common class of likelihood. However, in paleo-environmental studies, it is typical for proxy data to be subject to complex likelihoods (Ashe et al., 2022; Hibbert et al., 2016). For instance, organic matter that has been radiocarbon-dated undergoes a calibration procedure to account for the time-evolving atmospheric carbon concentration, which can yield a data chronology characterized by multi-modal distributions that significantly differ from each other. Similarly, it is common for paleo-environmental studies to use multiple types of proxy data with different likelihoods to infer a common signal. Recently, new approaches have been developed to account for non-parametric proxy distributions within a hierarchical modeling framework (e.g. Ashe et al., 2022), which could better characterize the underlying uncertainty but can be computationally expansive.

While PaleoSTeHM allows users to specify any number of change points (m in Eq. 10) in the model, determining the optimal number of change points can be challenging and may require additional modeling strategies. Recent advancements, such as Bayesian transdimensional models (e.g. Sambridge, 2016; Bodin et al., 2012; Gallagher et al., 2011), provide a flexible framework by treating the number of change points as an unknown parameter, allowing it to be inferred alongside other model parameters. Incorporating such approaches into PaleoSTeHM is a potential avenue for future development to address this complexity in abrupt paleo-environmental change problems.

The current GP kernel module incorporates commonly used kernel options that are stationary, isotropic, and space–time separable. While these assumptions simplify calculations significantly, they may not be suitable for some environmental applications. For example, temperature and dew point variations often exhibit strong non-stationary behavior influenced by diverse geographic and atmospheric conditions (Poppick and Stein, 2014). Additionally, the assumption of stationarity may cause rate uncertainty estimates to fail in properly reflecting the reduced uncertainty expected during periods with abundant data, as shown in Fig. 5d, largely because the model's variance is uniformly applied across the entire temporal domain (Heinonen et al., 2016). Furthermore, temperature anomalies over the last 2 millennia (Mann et al., 2008) demonstrate strong space–time interactions, which cannot be captured by a space–time-separable kernel (Tingley et al., 2012). Developing a scientifically richer class of kernel structures could be an important future advancement for PaleoSTeHM. However, given the fundamental differences across various paleo-environmental problems, generalizing sophisticated kernel structures to multiple fields remains challenging.

Another outstanding issue for GP-based process-level models is scalability: the standard GP models included in PaleoSTeHM v1.0 cannot scale well to large datasets (>10 000 data points) due to the computational cost, which increases at a rate of 𝒪(n³), where n is the number of data points (Hensman et al., 2013). Thus, implementing alternative classes of GP models within PaleoSTeHM to model large datasets, especially when incorporating modern environmental observations, which often consist of millions of data points, is an important next step for PaleoSTeHM to develop in the future. Some potentially efficient methods include sparse GP (Quinonero-Candela and Rasmussen, 2005), stochastic variational GP (Hensman et al., 2013), and exact GP with black-box matrix–matrix inference (Wang et al., 2019).

Building upon machine learning infrastructure, another promising direction for the future development of PaleoSTeHM is integrating spatiotemporal hierarchical modeling with machine-learning-based emulators as a process-level model. An emulator indicates a statistical model that mimics the behavior of the physics-based simulator but is computationally cheap to run (Reichstein et al., 2019), which is particularly useful for fast sensitivity analysis, model parameter calibration, and derivation of confidence intervals for the estimate. The use of statistical emulators trained by physical models will enable hierarchical models to capture the non-stationary physical systems better and enable better interpretation of the modeling results. For paleo-environment, Holden et al. (2019) present a GP-based emulator for an atmosphere–ocean general circulation model with intermediate complexity, and Lin et al. (2023b) developed a neural-network-based emulator for GIA-induced global sea-level change.

More broadly, PaleoSTeHM has been developed by a small team specialized in modeling paleo-sea-level changes over multi-millennial timescales. Moving forward, a critical objective is to expand PaleoSTeHM into a larger-scale paleo-environmental community project, where modules are developed autonomously by diverse research teams. The design of PaleoSTeHM, which allows modules to act as wrappers for independently developed code, is specifically intended to facilitate this collaborative effort.

6 Conclusions

Paleo-environmental records provide critical out-of-sample information essential for contextualizing current global changes and testing models used to simulate future environmental scenarios. However, our understanding of past environmental changes is often complicated by the sparse nature of geological records, geochronological uncertainties, and the indirect relationships between proxies and ecological variables. Hierarchical modeling offers a conceptually straightforward framework to address these challenges, though the limited availability of user-friendly software often hinders it. PaleoSTeHM offers a flexible and open-source platform that facilitates the rapid and easy implementation of hierarchical models for paleo-environmental applications. The inclusion of multiple process-level models in PaleoSTeHM allows it to be readily applicable across a broad spectrum of paleo-environmental studies. Additionally, its flexibility allows customization to meet the specific needs of diverse paleo-environmental problems, such as using different Gaussian process kernels or substituting alternative process-level models.

Appendix A: Additional model information

Table A1 Summary of model characteristics. The posterior is reported with a mean value with 90 % credible interval. GBR: Great Barrier Reef. NJ: New Jersey, NC: North Carolina. EIV: errors-in-variables. NI: noisy input.

Figure A1 Prior and posterior predictive checks for modeling results presented in the main text, Sect. 4.1.1. A random data point is selected for illustrative purposes. (a) Prior predictions compared with observational data, assuming uniform (blue) and normal (orange) likelihood functions. (b) Posterior predictions compared with observational data, assuming a uniform likelihood function. (c) Posterior predictions compared with observational data, assuming a normal likelihood function.

Figure A2 Residual plots for the three process-level models introduced in main text, Sect. 4.1.2. The weighted mean squared error (wMSE) for each model is given above each plot. Dashed red lines represents 0 error.

Figure A3 Optimization trace plot of empirical Bayesian analysis introduced in main text, Sect. 4.1.3.

Figure A4 Optimization trace plot of variational Bayesian analysis introduced in main text, Sect. 4.1.3. Note that the parameters shown here are variational inference parameters in Pyro, which are optimized to approximate the posterior distribution but do not directly correspond to the actual parameters of probability distribution.

Figure A5 Model prediction from model (ii) for relative sea level at time point −5500 CE along the US Atlantic coast. (a) Mean relative sea-level prediction, representing the sum of components in panels (b) and (c). (b) Relative sea-level prediction derived from the Gaussian process mean function, based on a glacial isostatic adjustment model incorporating the ICE_7G ice model and VM5a Earth model (Roy and Peltier, 2018; Peltier et al., 2015). (c) Relative sea-level change induced by the GP covariance function. (d) Standard deviation of relative sea level, as estimated by the GP covariance function.