Western Greenland ice sheet retreat history reveals elevated precipitation during the Holocene thermal maximum

We investigate changing precipitation patterns in the Kangerlussuaq region of western central Greenland during the Holocene thermal maximum (HTM), using a new chronology of ice sheet terminus position through the Holocene and a novel inverse modeling approach based on the unscented transform (UT). The UT is applied to estimate changes in annual precipitation in order to reduce the misfit between modeled and observed terminus positions. We demonstrate the effectiveness of the UT for time-dependent data assimilation, highlighting its low computational cost and trivial parallel implementation. Our results indicate that Holocene warming coincided with elevated precipitation, without which modeled retreat in the Kangerlussuaq region is more rapid than suggested by observations. Less conclusive is whether high temperatures during the HTM were specifically associated with a transient increase in precipitation, as the results depend on the assumed temperature history. Our results highlight the important role that changing precipitation patterns had in controlling ice sheet extent during the Holocene.

1 Introduction

During the early Holocene (∼11.7–8 ka), terrestrial and marine climate proxies from the Northern Hemisphere reveal a warmer than present peak in temperature (Kaufman et al., 2004; Marcott et al., 2013). This period of elevated temperatures, likely initiated by greater than modern insolation, is referred to as the Holocene thermal maximum (HTM). Its onset, duration, and severity were likely spatially variable (Kaufman et al., 2004). Records of HTM warming can be found in Greenland ice core records. For example, temperatures measured in the Dye-3 borehole show a pronounced HTM signal occurring from 7 to 4 ka and having values 2.5^∘C above present temperatures (Dahl-Jensen et al., 1998; Miller et al., 2010), whereas at the GISP2 site, the HTM appears to occur slightly earlier, following the 8.2 ka cold event (Kobashi et al., 2017) (Fig. 1).

While warming during the HTM is well established, less is known about the regional changes in precipitation that accompanied increased temperatures. Ice core records provide long-term estimates of accumulation (Alley et al., 1993), but these point measurements near ice divides are not representative of the precipitation across the ice sheet, particularly at lower elevations near the coast. Because the HTM was accompanied by lower Arctic sea ice extent (Polyak et al., 2010), it is possible that additional moisture was available to the Greenland ice sheet (GrIS) from open Arctic waters. This is supported by proxy evidence showing an increase in winter precipitation in western Greenland coincident with HTM warming (Thomas et al., 2016). However, temperature is known with greater certainty than precipitation.

Understanding feedbacks between temperature and precipitation during the HTM has implications for the future of the GrIS. Warming and declining sea ice are projected to cause an increase in Arctic precipitation (Bintanja and Selten, 2014; Singarayer et al., 2006). On a global scale, the moisture content of the atmosphere increases by around 7 % for every degree of warming, according to the Clausius–Clapeyron relation. On a regional scale, declining Arctic sea ice is expected to cause changes in atmospheric circulation, bringing more moisture to the Arctic (Bintanja and Selten, 2014). While there are important differences between HTM and modern climate, the history of retreat in western Greenland may provide insights into how the GrIS will respond to a warmer and possibly wetter future climate.

Modeling studies indicate that Holocene retreat in land-terminating regions of the GrIS were controlled primarily by surface mass balance rather than ice dynamics (Cuzzone et al., 2019; Lecavalier et al., 2014). Given the primary importance of surface mass balance in controlling modeled retreat, we explore the hypothesis that enhanced winter snowfall during the HTM may have slowed retreat by partially offsetting increased surface melt (Thomas et al., 2016). We investigate changes in precipitation in a land-terminating sector of the western central GrIS, near Kangerlussuaq, taking advantage of a new chronology of ice sheet terminus position (Young et al., 2020) and a novel inverse modeling approach based on the unscented transform (UT) (Julier and Uhlmann, 1997). In particular, we use the UT to estimate changes in annual precipitation during the Holocene by reducing the misfit between modeled and observed terminus positions in an isothermal flow line ice dynamics model (Brinkerhoff et al., 2017) (Sect. 2.1).

The inverse problem is posed as a Bayesian inference problem, and its solution involves estimating a non-Gaussian posterior probability distribution. Markov chain Monte Carlo (MCMC) methods, such as the Metropolis–Hastings method (Chib and Greenberg, 1995), provide one means of solving the inference problem by generating random samples from the posterior distribution. Generating samples from the posterior, however, requires repeatedly running the ice dynamics model with different precipitation histories as input, which is intractable even for a relatively computationally inexpensive flow line model.

The unscented transform provides a computationally efficient and trivially parallelizable alternative to MCMC methods. The basic idea of the unscented transform is to use a small, fixed number of deterministic sample points in order to estimate the statistical moments (e.g., mean and covariance) of the posterior distribution. Sigma points, each of which represents a different precipitation history input to the ice dynamics model, are generated a priori (Sect. 2.5.3). Consequently, all model runs can be performed simultaneously in parallel resulting in at least a 100-fold speed-up compared to MCMC methods.

2 Numerical methods for inference

2.1 Ice sheet model

We use the 1D isothermal flow line model with higher-order momentum balance described in Brinkerhoff et al. (2017). The momentum conservation equations are simplified using the Blatter–Pattyn approximation, assuming hydrostatic pressure and negligible vertical resistive stresses (Blatter, 1995; Pattyn, 2003). Default parameter values used in this work are specified in Table 1.

We adopt a linear sliding law of the form

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{b}}={\mathit{\beta }}^{\mathrm{2}}N{u}_{\mathrm{b}},\end{array}$$

where τ_b is basal shear stress, β² is a constant basal traction parameter, N is effective pressure, and u_b is sliding speed. Based on borehole water pressure measurements in Wright et al. (2016), basal water pressure P_w is assumed to be a fixed fraction P_frac=0.85 of ice overburden pressure P₀. Effective pressure is therefore given by

$$\begin{array}{}\text{(2)}& N=P_{\mathrm{0}}-P_{\text{frac}}{P}_{\mathrm{0}}.\end{array}$$

The basal traction parameter β² is tuned to minimize the misfit between modeled and observed surface velocities from Mouginot et al. (2017) for present-day Isunnguata Sermia values.

2.2 Flow line selection and moraine age constraints

To define the path followed by ice, we assume that flow follows the modern surface velocity field inland of the present-day margin. In ice-free regions, the direction of ice flow is inferred from bedrock topography (Fig. 1). Since the direction of ice flow is unknown and time varying, we cannot directly quantify the uncertainty introduced by errors in flow line selection. To account for some of this uncertainty, we perform inversions on two plausible adjacent paleo-flow lines in the Kangerlussuaq area.

The rate of Holocene retreat on each flow line is estimated using constraints on ice sheet terminus position from (Young et al., 2020) (Fig. 1). Terminus position data in Young et al. (2020) indicate that in the early Holocene (11.6 ka), the ice sheet margin was tens of kilometers inland of the present-day coastline. Although the moraine patterns are spatially complex, generally speaking there was a period of moderate retreat (∼10 km on the northern flow line and ∼30 km on the southern flow line) from 11.6 to 10.3 ka, followed by rapid retreat (∼100 km on both flow lines) from 10.3 to 8.1 ka. By 8.1 ka, the margin position was within 20 km of its present position on both flow lines (Fig. 5). The modern terminus position provides one additional constraint. Moraine ages have uncertainties of up to ±400 years.

For modern bedrock geometry along the flow lines, we use BedMachine v3 (Morlighem et al., 2017). Isostatic uplift and relative sea level changes are accounted for using a glacial isostatic adjustment model (Caron et al., 2018). This model, combined with the retreat chronology in Young et al. (2020), indicates that ice remained grounded on both the northern and southern flow lines from 11.6 ka onward.

2.3 Positive degree day model

Surface mass balance is estimated using a positive degree day (PDD) model (Johannesson et al., 1995). Annual surface mass balance is constructed in the PDD model using estimates of average monthly precipitation and temperature. Inputs into the PDD model include the unknown ice surface elevation S, modern monthly temperature T_m and precipitation P_m along the flow lines, as well as the seasonal temperature anomaly ΔT. As in Cuzzone et al. (2019), modern temperature and precipitation are computed as 30-year averages from 1980 to 2010 using the method of Box (2013).

To assess the sensitivity of modeled retreat history to temperature, we perform experiments using temperature reconstructions from both Buizert et al. (2018) and Dahl-Jensen et al. (1998). For the spatially explicit Buizert et al. (2018) reconstruction, monthly temperature anomalies are computed as averages along the flow lines. In contrast, Dahl-Jensen et al. (1998) reconstruct temperature only at the GRIP and Dye-3 borehole locations (Fig. 1). Since a full Holocene reconstruction is unavailable at the Dye-3 borehole site, which is closer to Kangerlussuaq, we use the temperature reconstruction at GRIP. The HTM is roughly 0.3 ^∘C warmer and 1500 years later in Dahl-Jensen et al. (1998) than in Buizert et al. (2018) (Fig. 6).

A limitation of the Dahl-Jensen et al. (1998) reconstruction is that it does not resolve seasonal temperatures. To address this, we calculate the difference between monthly and mean annual temperatures in Buizert et al. (2018) and apply those offsets to the mean annual temperature at GRIP from Dahl-Jensen et al. (1998).

Surface temperature T is computed monthly as

$$\begin{array}{}\text{(3)}& T=T_{\mathrm{m}}+\mathrm{\Delta }T+\mathit{\alpha }(S-S_{\mathrm{m}}),\end{array}$$

where S_m is the modern surface elevation and α= 5^∘C km⁻¹ is the lapse rate (Abe-Ouchi et al., 2007). Following Ritz et al. (2001) and Cuzzone et al. (2019), precipitation P along the flow line is estimated based on the Clausius–Clapeyron relation. In particular, precipitation is estimated by

$$\begin{array}{}\text{(4)}& P=P_T+\mathrm{\Delta }P=P_{\mathrm{m}}\exp \left({\mathit{\lambda }}_P(T-T_{\mathrm{m}})\right)+\mathrm{\Delta }P.\end{array}$$

The term P_T accounts for changes in precipitation solely due to changes in temperature. Here λ_P=0.07, which results in a 7 % increase in precipitation for every 1 ^∘C increase in temperature above modern (Abe-Ouchi et al., 2007; Ritz et al., 2001).

The term P_T does not capture the effects of many unknown climate factors that may have caused dynamic regional changes in Holocene precipitation. Therefore, we introduce a precipitation anomaly term ΔP, analogous to the temperature anomaly ΔT. This time-dependent function, which has units of meters water equivalent per annum (m w.e. a⁻¹), is used to adjust precipitation uniformly across a flow line in order to reduce mismatch between modeled and observed terminus positions. Unlike ΔT, which can be inferred from ice cores, ΔP will be used as a control variable to be determined using the inverse methods detailed in Sect. 2.5.3. Equations (3) and (4) provide a method of accounting for elevation changes through time and downscaling inputs to match the mesh resolution of the model (∼1 km).

Positive degree days and snowfall are computed month by month based on mean monthly temperature and precipitation (Johannesson et al., 1995). Snow is melted first at a rate of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ m w.e. ^∘ d⁻¹ followed by ice at a rate of $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{3}}$ m w.e. ^∘ d⁻¹. Snow melt is initially supposed to refreeze in the snowpack as superimposed ice. Runoff begins when the superimposed ice reaches a given fraction (60 %) of the snow cover (Reeh, 1991). A listing of ice flow and PDD model parameters is provided in Table 1, and all data sets used in the model are shown in Table 2.

2.4 Modeling limitations

A limitation of our modeling approach is that we do not account for potential ice dynamical effects caused by changes in surface runoff or subglacial hydrology. Modeling melt water runoff would be difficult in a flow line model due to flux of melt water in and out of the path of ice flow. Another limitation is that our model is isothermal. Unless ice temperature is treated in a vertically averaged sense, resolving temperature would require a 2D mesh, which would considerably increase the computational cost of the model. We consider the consequences of this simplification in Sect. 3.3, where we test sensitivity to the ice hardness parameter.

The PDD scheme outlined in Sect. 2.3 does not account for changes in surface mass balance due to orographic forcing or other complex interactions between the ice sheet and climate system that could be captured by coupling the ice sheet model to an Earth system model (Bahadory and Tarasov, 2018). Since our emphasis is on estimating precipitation and surface mass balance using an inverse modeling approach, we believe the computational cost of such an approach outweighs the benefits. Uncertainties related to feedbacks between ice dynamics and climate are assessed via extensive sensitivity testing (Sect. 3.3).

Figure 1 K marks the position of Kangerlussuaq, while D and G mark the locations of the Dye-3 and GISP2 boreholes, respectively. Northern and southern paleo-flow lines are shown as blue and red lines running left to right. The inset shows a detailed view of the modeled region. Modern bedrock elevation is expressed in meters above sea level. Historical moraines dating from 11.6 to 7.2 ka are indicated by labeled lines.

Table 1 Summary of primary model parameters used in this work. Default values are provided where applicable.

2.5 Data assimilation approach

In order to assess the initial mismatch between modeled and observed retreat histories, we perform a reference experiment with ΔP=0 and ΔT estimates from both Buizert et al. (2018) and Dahl-Jensen et al. (1998). To improve the fit to observations, we assimilate terminus position data to obtain improved estimates of Holocene precipitation anomalies. Previous modeling studies indicate that Holocene retreat in land-terminating sectors of the GrIS were dominated by surface mass balance rather than ice dynamics (Cuzzone et al., 2019; Lecavalier et al., 2014). Uncertainty in Holocene climate, and consequently surface mass balance, is therefore likely the primary cause of discrepancies between modeled and observed terminus positions.

In principle, ΔP, ΔT, or both could be tuned to improve the fit between modeled and observed terminus positions. We focus on precipitation because it is more poorly constrained than temperature. In the upcoming sections, we introduce a framework for time-dependent data assimilation based on the unscented transform (UT). Sections 2.5.1–2.5.2 outline the basic tenets of the UT. Sections 2.5.3–2.5.7 outline how the UT can be applied to estimate precipitation anomalies.

2.5.1 Overview of the unscented transform

In what follows, the notation x∼𝒩(x₀, P_x) means that x is a normally distributed random variable with mean vector x₀ and covariance matrix P_x. Suppose that $\mathit{x}\sim \mathcal{N}({\mathit{x}}_{\mathrm{0}},{\mathbf{P}}_x$), and $\mathcal{F}:{\mathbb{R}}^{\mathrm{n}}\to {\mathbb{R}}^{\mathrm{m}}$ is a nonlinear function. We would like to estimate the distribution of the non-Gaussian random variable

$$\begin{array}{}\text{(5)}& \mathit{y}=\mathcal{F}\left(\mathit{x}\right)+\mathit{ϵ},\end{array}$$

where $\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathbf{R}$) is the measurement noise.

In general, the non-Gaussian probability distribution for y can be approximated using Markov chain Monte Carlo (MCMC) methods such as the Metropolis–Hastings algorithm (Chib and Greenberg, 1995). However, if the nonlinear function is time-consuming to compute, generating thousands of MCMC samples is often intractable. As a computationally efficient alternative to MCMC methods, Julier and Uhlmann (1997) introduced a method for approximating the mean and covariance of y called the unscented transform (UT) ¹.

The term unscented transform has been applied somewhat broadly to a family of methods that approximate the statistical moments of a non-Gaussian random variable using a small, deterministic set of sample points called sigma points. It is known primarily in the context of the unscented Kalman filter. However, the UT can be applied more generally as an alternative to traditional MCMC methods. Sigma points and weight sets are designed to accurately estimate moments of a transformed random variable using a minimal number of function evaluations.

A set of vectors, called sigma points, are chosen with the same weighted sample mean and weighted covariance structure as x. There are many algorithms for generating sigma points sets with different numbers of points and orders of accuracy. A commonly used set of 2n+1 sigma points is given by

$$\begin{array}{}\text{(6)}& {\mathit{\chi }}_i=\left\{\begin{array}{ll}{\mathit{x}}_{\mathrm{0}}& i=\mathrm{0}\\ {\mathit{x}}_{\mathrm{0}}-\sqrt{n+\mathit{\kappa }}[\sqrt{{\mathbf{P}}_x}{]}_i& i=\mathrm{1},\mathrm{\dots },n\\ {\mathit{x}}_{\mathrm{0}}+\sqrt{n+\mathit{\kappa }}[\sqrt{{\mathbf{P}}_x}{]}_i& i=n+\mathrm{1},\mathrm{\dots },\mathrm{2}n\end{array}\right.,\end{array}$$

with corresponding mean and covariance weights given by

$$\begin{array}{}\text{(7)}& w_i^{\left(\mathrm{m}\right)}=w_i^{\left(\mathrm{c}\right)}=\left\{\begin{array}{ll}\mathit{\kappa }/(n+\mathit{\kappa })& i=\mathrm{0}\\ \mathrm{1}/\mathrm{2}(n+\mathit{\kappa })& \text{otherwise}\end{array}\right..\end{array}$$

The notation ${\left[\sqrt{{\mathbf{P}}_x}\right]}_i$ refers to the ith row of a matrix square root (typically computed by Cholesky factorization) of P_x, and κ is a free parameter controlling the scaling of the sigma points around the mean. Julier and Uhlmann (1997) recommend a default value of $\mathit{\kappa }=\mathrm{3}-n$. However, κ can be fine-tuned to reduce prediction errors for a given problem.

Young et al. (2020)Morlighem et al. (2017)Mouginot et al. (2017)Box (2013)Buizert et al. (2018)Dahl-Jensen et al. (1998)Caron et al. (2018)

Table 2 Citations for the primary data sets used in this work.

The nonlinear function ℱ is applied to each sigma point to yield a set of transformed points

$$\begin{array}{}\text{(8)}& {\mathcal{Y}}_{\mathit{i}}=\mathcal{F}\left({\mathit{\chi }}_i\right).\end{array}$$

The mean $\stackrel{\mathrm{‾}}{\mathit{y}}$ and covariance matrix P_y of y are then estimated as weighted sums.

$$\begin{array}{}\text{(9)}& \stackrel{\mathrm{‾}}{\mathit{y}}=\sum _{i=\mathrm{0}}^{\mathrm{2}n}{w}_i^{\left(\mathrm{m}\right)}{\mathcal{Y}}_{\mathit{i}}\end{array}$$

$$\begin{array}{}\text{(10)}& {\mathbf{P}}_y=\sum _{i=\mathrm{0}}^{\mathrm{2}n}{w}_i^{\left(\mathrm{c}\right)}({\mathcal{Y}}_{\mathit{i}}-\stackrel{\mathrm{‾}}{\mathit{y}})({\mathcal{Y}}_{\mathit{i}}-\stackrel{\mathrm{‾}}{\mathit{y}}{)}^T+\mathbf{R}\end{array}$$

A visual example of this algorithm is shown in Fig. 2.

Figure 2 (a) White points are samples from a 2D Gaussian distribution $\mathit{x}\sim \mathcal{N}(\mathbf{0},I$). The red ellipse represents the covariance of the distribution, while red points are sigma points χ_i used for the unscented transform (UT). (b) White points are samples from the transformed, non-Gaussian distribution for $\mathit{y}=\mathcal{F}\left(\mathit{x}\right)+\mathit{ϵ}$ with $\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathbf{R}$). Red points are transformed sigma points 𝓨_i=ℱ(χ_i). UT approximations of the mean $\stackrel{\mathrm{‾}}{\mathit{y}}$ and covariance P_y of y are estimated via weighted sums of transformed sigma points.

2.5.2 Bayesian inference using the UT

Given a measurement of y_o, we would like to estimate the posterior distribution

$$\begin{array}{}\text{(11)}& \mathbf{P}\left(\mathit{x}\right|{\mathit{y}}_{\mathrm{o}})\propto \mathbf{P}(\mathit{y}\left|\mathit{x}\right)\mathbf{P}\left(\mathit{x}\right).\end{array}$$

Using an approach called statistical linearization (Särkkä, 2013) the joint distribution for [x, y]^T can be approximated by

$$\begin{array}{}\text{(12)}& \left[\begin{array}{c}\mathit{x}\\ \mathit{y}\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}\mathit{x}\\ \mathit{\mu }\end{array}\right],\left[\begin{array}{cc}{\mathbf{P}}_x& {\mathbf{P}}_{xy}\\ {\mathbf{P}}_{xy}^T& {\mathbf{P}}_y\end{array}\right]\right),\end{array}$$

with

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\mu }=\sum _i{w}_i^{\left(\mathrm{m}\right)}{\mathcal{Y}}_{\mathit{i}},\text{(14)}& {\displaystyle }{\mathbf{P}}_y=\sum _i{w}_i^{\left(\mathrm{c}\right)}({\mathcal{Y}}_{\mathit{i}}-\mathit{\mu })({\mathcal{Y}}_{\mathit{i}}-\mathit{\mu }{)}^T+\mathbf{R},\text{(15)}& {\displaystyle }{\mathbf{P}}_{xy}=\sum _i{w}_i^{\left(\mathrm{c}\right)}({\mathcal{X}}_{\mathit{i}}-{\mathit{x}}_{\mathrm{0}})({\mathcal{Y}}_{\mathit{i}}-\mathit{\mu }{)}^T.\end{array}$$

Here, 𝓧_i and 𝓨_i=ℱ(𝓧_i) are sigma points and transformed sigma points, respectively. Matrices P_y and P_xy are known as the measurement covariance and cross covariance, respectively.

Given a measurement y₀, the joint distribution then easily yields a Gaussian approximation of the posterior distribution. Using the following equation:

$$\begin{array}{}\text{(16)}& \mathbf{K}={\mathbf{P}}_{xy}{\mathbf{P}}_y^{-\mathrm{1}},\end{array}$$

we have

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathit{x}}^{\mathbf{\prime }}={\mathit{x}}_{\mathbf{o}}+\mathbf{K}\left[{\mathit{y}}_{\mathrm{o}}-\mathit{\mu }\right],\text{(18)}& {\displaystyle }{\mathbf{P}}^{\prime }={\mathbf{P}}_x-{\mathbf{KP}}_y{\mathbf{K}}^T,\end{array}$$

where x^′ and P^′ are approximations of the posterior mean and covariance, respectively. Readers familiar with Kalman filters might recognize that x^′ and P^′ are computed using a Kalman “update” step given a measurement y_o and Kalman gain K (Särkkä, 2013).

2.5.3 Assimilating glacier length observations

Time-dependent data assimilation using the UT involves running the ice sheet model with a set of different precipitation anomaly histories, each corresponding to a different sigma point. This is followed by a post-processing step, which incorporates the ice sheet terminus chronology data via a correction of the prior mean vector and covariance matrix. Implementation of the unscented transform is straightforward and easily parallelizable since each model run is independent. In the following section, we outline the mathematical details of this process.

We seek to find ΔP histories that match the observed retreat history on both flow lines. An optimal solution should reproduce the observed retreat history within uncertainty while not overfitting the data. We discretize the problem by estimating the precipitation anomaly at times $t_{\mathrm{1}},t_{\mathrm{2}},\mathrm{\dots },t_n$ ranging from 11.6 to 0 ka. In practice, we use a regular grid of 44 points, spaced roughly 250 years apart. Precipitation anomaly values at these time points are denoted by $\mathrm{\Delta }{p}_{\mathrm{1}},\mathrm{\Delta }{p}_{\mathrm{2}},\mathrm{\dots },\mathrm{\Delta }{p}_n$, respectively, and assembled in a vector

$$\begin{array}{}\text{(19)}& \mathbf{\Delta }\mathit{p}=[\mathrm{\Delta }{p}_{\mathrm{1}},\mathrm{\Delta }{p}_{\mathrm{2}},\mathrm{\dots },\mathrm{\Delta }{p}_n{]}^T.\end{array}$$

Given a multivariate Gaussian prior $\mathbf{\Delta }\mathit{p}\sim \mathcal{N}(\mathbf{\Delta }{\mathit{p}}_{\mathrm{0}},{\mathbf{P}}_{\mathrm{0}}$), which encodes assumptions about the structure of the precipitation anomaly (Sect. 2.5.4), we would like to estimate the mean and covariance of the posterior distribution

$$\begin{array}{}\text{(20)}& P\left(\mathbf{\Delta }\mathit{p}\right|{\mathbf{\ell }}_{\mathrm{0}})\propto P({\mathbf{\ell }}_{\mathrm{0}}\left|\mathbf{\Delta }\mathit{p}\right)P\left(\mathbf{\Delta }\mathit{p}\right).\end{array}$$

Here, the measurement vector

$$\begin{array}{}\text{(21)}& {\mathbf{\ell }}_{\mathrm{0}}=[{\mathrm{\ell }}_{\mathrm{1}},{\mathrm{\ell }}_{\mathrm{2}},\mathrm{\dots },{\mathrm{\ell }}_{\mathrm{m}}{]}^T,\end{array}$$

contains measured glacier lengths at discrete points in time. Our procedure for defining the measurement mean y_o and covariance R are discussed in Sect. 2.5.5.

We can think of the ice sheet model as a function that maps precipitation anomaly inputs to glacier length outputs with some additive observation noise

$$\begin{array}{}\text{(22)}& \mathbf{\ell }=\mathcal{F}\left(\mathbf{\Delta }\mathit{p}\right)+\mathit{ϵ}.\end{array}$$

Discrete precipitation anomaly values are linearly interpolated for input into the ice dynamics model, which has time steps on the order of months. The function ℱ returns glacier lengths at the same m discrete times as in ℓ₀.

To predict the posterior distribution, we use the methods outlined in Sect. 2.5.1 and 2.5.2. First, sigma points are generated based on the prior distribution for Δp. To reduce computational costs, we use a minimal set of $n+\mathrm{1}=\mathrm{45}$ sigma points 𝓟_i with corresponding weights $w_i^{\left(\mathrm{m}\right)}=w_i^{\left(\mathrm{c}\right)}$ generated using the method presented in Menegaz et al. (2011). Their method includes one free parameter $\mathrm{0}<w_{\mathrm{0}}<\mathrm{1}$, which can be tuned to reduce prediction errors. While this method has a lower order of accuracy than other methods, we find it often produces comparable results to other larger sigma point sets in practice.

Sigma points are propagated through the model to obtain transformed points 𝓛_i=ℱ(𝓟_i). In this context, sigma points 𝓟_i correspond to different time-dependent precipitation anomaly histories, while the transformed points 𝓛_i correspond to the resulting glacier length histories when those precipitation anomalies are used as inputs (Fig. 3). The structure of the sigma points reflects the mathematical formulation of the Menegaz et al. (2011) sigma points. Hence, they are not merely random samples from the prior distribution.

Transformed sigma points are computed simultaneously in parallel, using one core per sigma point. After all transformed sigma points have been computed, the mean and covariance of the posterior distribution are estimated as outlined in Sect. 2.5.2. In parallel, this procedure takes roughly the same amount of time as a single forward model run. Unlike a standard filtering approach to data assimilation, all measurement data are incorporated simultaneously rather than time step by time step. For that reason, the Kalman update step corrects the entire time-dependent precipitation history at once. Moreover, unlike in Kalman smoothing, we approximate the full posterior distribution rather than the probability distributions

$$\begin{array}{}\text{(23)}& P\left(\mathrm{\Delta }{p}_i\right|{\mathbf{\ell }}_{\mathrm{0}})i=\mathrm{1},\mathrm{\dots },n.\end{array}$$

Note that the variables Δp_k with k≠i are marginalized out of the Kalman smoothing distributions. The use of time-dependent sigma points distinguishes our approach from standard Kalman filtering or Kalman smoothing approaches, and does not rely on the assumptions that states (Δp_i) and measurements (ℓ_i) satisfy the Markov property.

2.5.4 Gaussian process prior for regularization

We adopt a Gaussian process prior (Rasmussen, 2004) to control the temporal smoothness of ΔP. A Gaussian process can be thought of as a distribution over functions. That is, random samples from a Gaussian process are functions rather than individual points or vectors. A collection of random variables $\mathit{\left\{}f\right(t):t\in \mathcal{T}\mathit{\}}$ is said to be drawn from a Gaussian process with mean function m(⋅) and covariance function $k(\cdot ,\cdot$) if, for any finite set of elements $t_{\mathrm{1}},\mathrm{\dots },t_n\in \mathcal{T}$, the random variables $f\left(t_{\mathrm{1}}\right),\mathrm{\dots },f(t_n$) have the distribution

$$\begin{array}{}\text{(24)}& \mathit{f}\sim \mathcal{N}\left(\mathit{m},\mathbf{K}\right),\end{array}$$

with

$$\begin{array}{}\text{(25)}& \mathit{f}=\left[f\right(t_{\mathrm{1}}),f(t_{\mathrm{2}}),\mathrm{\dots },f(t_n){]}^T,\end{array}$$

$$\begin{array}{}\text{(26)}& \mathit{m}=\left[m\right(t_{\mathrm{1}}),m(t_{\mathrm{2}}),\mathrm{\dots },m(t_n){]}^T,\end{array}$$

and

$$\begin{array}{}\text{(27)}& \mathbf{K}=\left[\begin{array}{ccc}k(t_{\mathrm{1}},t_{\mathrm{1}})& \mathrm{\dots }& k(t_{\mathrm{1}},t_n)\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ k(t_n,t_{\mathrm{1}})& \mathrm{\dots }& k(t_n,t_n)\end{array}\right].\end{array}$$

The set 𝒯 is called the index set, and specifies the domain of the Gaussian process. Here, the index set represents points in time.

The prior distribution for Δp has mean vector Δp₀ and covariance matrix P₀=K of the form shown in Eq. (27). We use a squared exponential covariance function

$$\begin{array}{}\text{(28)}& k(t,t^{\prime })={\mathit{\sigma }}_P^{\mathrm{2}}\exp \left({\displaystyle \frac{(t-t^{\prime }{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\tau }}^{\mathrm{2}}}}\right),\end{array}$$

where σ² is a scaling constant and τ is a characteristic timescale. Variables Δp_i and Δp_k are more highly correlated the closer they are in time. In effect, this acts as a form of temporal regularization, in which smooth precipitation history functions are preferred over less smooth ones. We discuss the choice of the mean Δp₀ in Sect. 3.2.

Figure 3 (a) A subset of 6 out of a total of 45 Menegaz et al. (2011) sigma points 𝓟_i representing different precipitation anomaly histories. (b) Six corresponding glacier length histories 𝓛_i or transformed sigma points obtained by inputting the sigma points 𝓟_i into the ice dynamics model ℱ.

2.5.5 Measurement mean and variance

Observations of terminus position are available roughly every 1000 years between 11.6 and 7.2 ka, with a gap from 7.2 ka to present. In contrast, model time steps are on the order of months. Due to these disparate timescales, we use the following procedure to estimate the measurement mean ℓ₀ and covariance matrix R on a timescale more appropriate for the ice sheet model.

We define a Gaussian process prior of “candidate” glacier length histories ℓ_P(t) as follows. The mean function $\widehat{{\mathrm{\ell }}_P}(t$) is obtained by linearly interpolating between glacier length observations. The Brownian covariance kernel for the Gaussian process is defined by

$$\begin{array}{}\text{(29)}& k_P(t,t^{\prime })={\mathit{\sigma }}_{\mathrm{\ell }}^{\mathrm{2}}min(t,t^{\prime }).\end{array}$$

Candidate retreat histories are generated by drawing random samples from the Gaussian process. The Brownian covariance kernel results in random glacier length histories that are somewhat noisy but correlated over shorter timescales (Fig. 4). Candidate length histories are resampled so that the mean moraine formation times and uncertainties match the observations in Young et al. (2020). Hence, highly implausible candidate histories are rejected. The average length and variance of this culled set of samples is computed at a series of time slices to obtain a plausible measurement mean ℓ₀ and diagonal measurement covariance matrix R (Fig. 4).

2.5.6 Iterative optimization procedure

Optimizations are conducted in multiple passes. In the first pass, the measurement covariance matrix R is multiplied by a factor of 0.25 so that the measurements are initially weighted more than the prior. This produces a reasonable fit to the data, even given a poor initial estimate of ΔP. The optimal precipitation anomaly from a given iteration is used as the prior mean in the next iteration. We use the same prior covariance matrix P₀ for regularization in each iteration. After two to three iterations, modeled and observed terminus positions match within measurement uncertainty (Sect. 3). In our experience, the results of iteration are not dependent on the choice of prior mean in the first iteration, but we find that convergence can be improved by choosing a sensible initial guess, as in Sect. 3.2.

Figure 4 The measurement mean ℓ₀ (solid black line) and 95 % confidence bands (gray shaded region) for the northern flow line are estimated by generating random retreat histories with the same mean moraine formation ages and variances as the observations. The green, blue, and orange lines represent four random plausible retreat histories. Red dots denote the estimated mean moraine formation ages, while red lines show 95 % confidence intervals.

2.5.7 Approach to sensitivity testing

As described, the data assimilation method accounts for measurement but not model uncertainty. It can easily be extended to account for uncertainties in the ice flow and PDD model parameters. We define an augmented state vector

$$\begin{array}{}\text{(30)}& \mathit{u}=[\mathbf{\Delta }\mathit{p},\mathit{\theta }{]}^T,\end{array}$$

where θ is a vector of scalar parameters including the natural logarithm of the rate factor for ice, the basal traction parameter, a parameter controlling precipitation scaling with temperature, and the PDD melt rate parameters for ice and snow

$$\begin{array}{}\text{(31)}& \mathit{\theta }=\left[\ln \right(A),{\mathit{\beta }}^{\mathrm{2}},{\mathit{\lambda }}_P,{\mathit{\lambda }}_i,{\mathit{\lambda }}_{\mathrm{s}}{]}^T.\end{array}$$

The prior distribution for the augmented state vector is given by

$$\begin{array}{}\text{(32)}& \mathit{u}\sim \mathcal{N}\left(\left[\begin{array}{c}\mathbf{\Delta }{\mathit{p}}_{\mathrm{o}}\\ {\mathit{\theta }}_{\mathrm{0}}\end{array}\right],\left[\begin{array}{cc}{\mathbf{P}}_{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& \mathrm{\Theta }\end{array}\right]\right),\end{array}$$

where θ₀ is the parameter mean vector and Θ is a diagonal matrix containing parameter variances.

The unscented transform is applied to the augmented function

$$\begin{array}{}\text{(33)}& \mathbf{\ell }=\widehat{\mathcal{F}}\left(\mathit{u}\right)+\mathit{ϵ},\end{array}$$

to obtain estimates of the joint distribution for $[\mathbf{\Delta }\mathit{p},\mathit{\theta },\mathbf{\ell }{]}^T$ and the conditional distribution for $[\mathbf{\Delta }\mathit{p},\mathit{\theta }{]}^T|{\mathbf{\ell }}_{\mathrm{0}}$. Since parameters are included as state variables, sigma points reflect a variety of precipitation histories and parameter sets. A model run for a particular sigma point is initialized from an appropriate steady state using the parameter set for that point.

2.6 Model initialization

Model runs are initialized by tuning the precipitation anomaly to obtain a steady state at 12.6 ka, with a margin position 5 km beyond the 11.6 ka moraine. We invert for a precipitation anomaly time series that forces a retreat of 5 km over 1000 years to obtain an initial ice sheet configuration with the correct terminus position at 11.6 ka (Fig. 5). This initialization procedure is intended to ease the ice sheet out of steady state in order to avoid strong transient effects at the beginning of model runs.

Figure 5 Initial ice sheet configurations for the northern (a) and southern (b) paleo-flow lines. Solid black and blue lines represent the bedrock elevation and ice surface, respectively. Moraine positions are indicated by arrows with associated ages expressed in thousands of years before present.

3 Results

3.1 Reference experiment

To assess the initial misfit between modeled and observed Holocene retreat, the model is forced with ΔT reconstructions from Buizert et al. (2018) and Dahl-Jensen et al. (1998) and a zero precipitation anomaly. Precipitation is scaled with temperature according to Eq. (4), neglecting possible influences from changing Arctic sea ice cover, atmospheric circulation, or other unknown climate factors. Modeled ice retreat is far more rapid than observed on both flow lines (Fig. 8). Colder temperatures during the early Holocene (Fig. 8a) lead to a somewhat more plausible retreat history using the Dahl-Jensen et al. (1998) temperature forcing versus the Buizert et al. (2018) forcing. However, by 8 ka, ice has retreated inland of the present-day margin in both reconstructions.

3.2 Precipitation anomaly inversions

We estimate precipitation anomalies on the northern and southern flow lines using both the Buizert et al. (2018) and Dahl-Jensen et al. (1998) temperature reconstructions. Given the rapid retreat in the reference experiment, we expect that a positive precipitation anomaly will be required to match observed terminus positions in the early Holocene. Therefore, in the first round of optimization, we assume a prior mean of the form

$$\begin{array}{}\text{(34)}& \mathrm{\Delta }P={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-\widehat{\mathit{\tau }}),\end{array}$$

where $\widehat{\mathit{\tau }}$ is a rescaled time variable that increases from zero at 11.6 ka to one at 0 ka. The results of the iterative optimization procedure are insensitive to the prior mean selected in the first iteration. Using a sigma point scaling parameter w₀=0.5 for the Menegaz et al. (2011) sigma point set ensures that a wide region around the mean is explored in each iteration.

Positive precipitation anomalies are predicted throughout most of the Holocene for both temperature reconstructions (Fig. 6). While differences between the northern and southern flow lines are relatively minor, there are significant differences in precipitation between the Buizert et al. (2018) and Dahl-Jensen et al. (1998) inversions. In the Buizert et al. (2018) inversion, the largest precipitation anomalies (up to 1 m w.e. a⁻¹) occur during the early Holocene. Precipitation remains relatively high during the HTM (10–6 ka) but dips before the 8.2 ka cold event. For the Dahl-Jensen et al. (1998) inversion, ΔP is relatively low during the early Holocene but increases during the HTM (8.5–3 ka). Unlike the Buizert et al. (2018) inversion, there is an evident trend between HTM warming and increased snowfall (Fig. 9).

Forcing the model with mean estimated precipitation anomalies yields plausible retreat histories on both flow lines. Modeled and observed retreat chronologies match within uncertainty (Fig. 8). A significant fraction of Holocene precipitation (typically >90 %) falls as snow. Hence, a positive precipitation anomaly can be interpreted directly as additional snowfall and snow accumulation. Average HTM snowfall is around 35 % higher than modern in both temperature reconstructions (Fig. 9). However, overall trends in Holocene snowfall differ between reconstructions.

Figure 6 (a) Buizert et al. (2018) and Dahl-Jensen et al. (1998) mean annual temperature anomalies are shown as blue and red lines, respectively. Gray and black lines, respectively, show average December, January, February (DJF) and June, July, August (JJA) temperature anomalies for the Buizert et al. (2018) reconstruction. Shaded blue and pink regions demarcate the extent of the HTM for the Buizert et al. (2018) and Dahl-Jensen et al. (1998) reconstructions. (b) Estimated precipitation anomaly histories for the Buizert et al. (2018) and Dahl-Jensen et al. (1998) inversions are shown by blue and red lines, respectively. Solid lines show estimates for the northern flow line, while dashed lines are for the southern flow line.

Figure 7 Average precipitation (both solid and liquid) on the northern flow line from the Buizert et al. (2018) reconstruction (black line), which is based on the TraCE-21ka coupled ocean–atmosphere general circulation model of the last deglaciation versus estimated precipitation derived in this work based on reconciling the modeled and observed retreat chronologies (blue line).

Figure 8 Modeled retreat histories for the northern (a) and southern (b) flow lines using a variety of ΔT and ΔP forcings. Blue and red lines show modeled retreat using Buizert et al. (2018) and Dahl-Jensen et al. (1998) temperature forcings, respectively. Dashed lines show the results of forcing the model with ΔP=0, while solid lines show the results of using optimal estimated precipitation anomalies. Black diamonds and bars show reconstructed terminus positions, along with 95 % confidence intervals. Shaded blue and pink regions demarcate the extent of the HTM for the Buizert et al. (2018) and Dahl-Jensen et al. (1998) reconstructions, respectively.

Figure 9 HTM snowfall averaged across the northern and southern flow lines for the Buizert et al. (2018) (solid blue line) and Dahl-Jensen et al. (1998) (solid red line) inversions. Shaded blue and pink regions demarcate the extent of the HTM for the Buizert et al. (2018) and Dahl-Jensen et al. (1998) reconstructions, respectively. Dashed lines show pre-HTM, HTM, and post-HTM average snowfall for each inversion.

3.3 Sensitivity testing

We assess the sensitivity of Holocene precipitation anomalies to modeling uncertainties by performing an HTM inversion using the methodology described in Sect. 2.5.7. To obtain accurate uncertainty estimates, we use a fifth-order accurate sigma point set based on Li et al. (2017) (Appendix A), as we find that the second-order Menegaz et al. (2011) set likely underestimates covariance. Inversions are performed on the northern flow line using the Buizert et al. (2018) temperature reconstruction. Model runs are initialized from steady states around at 10.5 ka, 500 years prior to the Buizert et al. (2018) HTM. Prior and posterior parameter values are reported in Table 3.

Table 3 Summary of primary model parameters used in this work. Default values are provided where applicable. Prior parameter uncertainties assumed for the sensitivity test are shown in the 2σ column.

Mean HTM precipitation anomalies are within 2 cm w.e. a⁻¹ of the inversion presented in Sect. 3.2 (Fig. 10). Parameter uncertainties contribute to uncertainty in ΔP. Overall, however, temperature uncertainty is far more significant than uncertainty in model parameters. Differences in model initialization do not significantly impact the results of the inversion. Although sensitivity tests are initialized from steady states at 10.5 ka, while inversions in Sect. 3.2 are initialized from a transient state at 11.6 ka, the mean HTM precipitation anomalies are nearly identical. Estimated posterior parameter values are also similar to the assumed prior values (Table 3).

Figure 10 The blue line and shaded region show the estimated mean ΔP and 95 % confidence bands, respectively, for the sensitivity test, accounting for uncertainty in a number of ice flow and PDD model parameters. The dashed black line shows the mean estimated ΔP from a previous inversion assuming no uncertainty in model parameters.

4 Discussion

We infer changes in Holocene precipitation in the Kangerlussuaq region of western Greenland using a new chronology of ice sheet terminus position from Young et al. (2020) and an inverse modeling procedure based on the unscented transform. We find that scaling precipitation with temperature via the Clausius–Clapeyron equation (Eq. 4) results in excessively fast retreat during the early Holocene for both the Buizert et al. (2018) and Dahl-Jensen et al. (1998) temperature reconstructions. Thus, temperature-driven changes in precipitation and accumulation alone are not sufficient to reproduce the observed pattern of retreat. Inversions show that adding precipitation throughout the Holocene and specifically during the HTM yields a good fit to observations (Figs. 6, 8). Since a large fraction of precipitation falls as snow, a positive precipitation anomaly can be interpreted fairly directly as an increase in accumulation beyond what would be expected from temperature changes alone.

There are considerable differences in predicted precipitation anomalies depending on the assumed temperature history. Inversions using the Buizert et al. (2018) temperature reconstruction show generally decreasing snowfall through the Holocene. Due to high snowfall in the early Holocene (nearly 100 % higher than modern) and a dip in HTM precipitation associated with the 8.2 ka cold event, average HTM snowfall roughly matches the overall Holocene average, which is about 35 % above modern. In contrast, inversions using the Dahl-Jensen et al. (1998) temperature reconstruction show a clear trend between HTM warming and increased snowfall (Fig. 9). This trend could be interpreted as a transient increase in snowfall due to reduced Arctic sea ice cover or changes in atmospheric circulation during the HTM.

A large positive ΔP correction during the early Holocene in the Buizert et al. (2018) inversion likely reflects the dependence of precipitation on temperature. Low temperatures during the early Holocene result in low precipitation. Without an additional moisture source to increase snowfall, the ice sheet retreats far more rapidly than observed. In the Dahl-Jensen et al. (1998) reconstruction, warming occurs more gradually in the early Holocene. Accumulation and ablation are more closely balanced, resulting in a smaller precipitation anomaly correction from 11.6 to 10 ka.

When interpreting results of precipitation inversions, it is important to consider that the Buizert et al. (2018) and Dahl-Jensen et al. (1998) reconstructions are obtained using different methodologies. The Greenland-wide Buizert et al. (2018) reconstruction is obtained by merging the TraCE-21ka coupled ocean–atmosphere general circulation model of the last deglaciation (Liu et al., 2009; He et al., 2013) with borehole temperature reconstructions at GISP2, NGRIP, and NEEM. Our estimated precipitation history in the Kangerlussuaq region is significantly different from TraCE-21ka, particularly during the early Holocene and parts of the HTM (Fig. 7). In contrast to Buizert et al. (2018), the Dahl-Jensen et al. (1998) temperature reconstruction at GRIP offers only a pointwise estimate containing no explicit spatial or seasonal information.

While the Buizert et al. (2018) temperature reconstruction is arguably more suitable for our purposes, since it resolves spatial and seasonal patterns in ΔT, there is still considerable uncertainty in temperature, particularly near the ice margins. Despite this, there is some consensus between the Buizert et al. (2018) and Dahl-Jensen et al. (1998) inversions. For example, average HTM snowfall is roughly 35 % higher than modern in both cases (Fig. 9).

Predicted HTM precipitation anomalies are not particularly sensitive to uncertainties in the ice sheet or PDD model parameters or the initialization procedure. In the sensitivity test presented in Sect. 3.3, model runs are initialized using steady states at 10.5 ka. Runs in the full Holocene inversions in Sect. 3.2, which use fixed parameter sets, are initialized from a transient state at 11.6 ka. Even considering these differences in model initialization, estimated HTM precipitation anomalies are comparable in all inversions (Fig. 10). Additionally, ΔP does not appear to be sensitive to uncertainties in bedrock geometry, as inversions conducted with a static bedrock geometry yielded similar results to inversions accounting for isostatic uplift.

5 The unscented transform as a data assimilation method

Significant strides have been made in time-dependent data assimilation in glaciology using adjoint-based methods. Goldberg and Heimbach (2013) infer the initial thickness and basal conditions for a synthetic ice sheet given snapshots of ice thickness at discrete times. Larour et al. (2014) demonstrate a data assimilation framework within the Ice Sheet System Model (ISSM), capable of obtaining temporal estimates of surface mass balance and basal friction given surface altimetry.

In contrast to adjoint-based approaches, the unscented transform (UT) does not require computing the Jacobian or Hessian of an objective function or a special checkpointing code for time-dependent problems. Adjoint-based methods are advantageous for extremely high-dimensional problems, as the required number of model runs is independent of the number of parameters. However, the UT provides more accurate uncertainty estimates than linearization (Julier and Uhlmann, 1997). Hessian information can be used to improve uncertainty estimates (Isaac et al., 2015). However, this methodology uses purely local approximations to the nonlinear function around the maximum a posteriori probability (MAP) estimator, which may affect the quality of uncertainty estimates.

The unscented transform has advantages over Markov chain Monte Carlo methods for inference problems with a relatively small number of unknown parameters (<1000 or so parameters). In this work, we optimize for n=44 parameters representing precipitation anomaly values at discrete points in time. Since function evaluations at each sigma point are independent, the UT is trivially parallelizable. Consequently, on a high-end desktop, the iterative optimization process presented in Sect. 2.5.6 takes roughly the same amount of time as performing three successive forward model runs.

MCMC methods can be parallelized to some extent by utilizing multiple interacting Markov chains running in parallel (e.g., Chowdhury and Jermaine, 2018) or by combining samples from independent chains in a post-processing step (e.g., Neiswanger et al., 2013). Nonetheless, computational bottlenecks persist in MCMC methods since individual Markov chains are inherently serial. Consequently, UT approximations to the posterior distribution can be generated hundreds to thousands of times faster for small parameter sets.

A drawback of the UT is that its accuracy is inherently limited by using a predetermined number of sample points. In our case, for example, using a minimal set of n+1 sigma points for n unknown parameters yields accurate mean estimates but appears to underestimate covariance when compared to a higher-order cubature method. Without testing multiple sigma point sets and scaling parameters, it can be difficult to assess the accuracy of UT estimates of the posterior. MCMC methods, in contrast, can provide arbitrarily accurate estimates of the posterior given sufficient computation time. Moreover, they can resolve the full posterior distribution rather than computing its moments as in the UT.

An alternative to traditional MCMC methods is to use surrogate models or emulators (Gong and Duan, 2017). Here, a computationally inexpensive surrogate model is trained to approximate the output of a more complex model function. The surrogate model can then be used in place of the full model for the purpose of MCMC sampling, significantly reducing the overall computational cost. A few types of surrogate models have already been applied to glaciological problems including deep neural networks (Tarasov et al., 2012) and Gaussian processes (e.g., Chang et al., 2016; Pollard et al., 2016).

Surrogate models require an initial training phase. Pollard et al. (2016), for example, perform 625 ice sheet model runs using combinations of four unknown simulation parameters. A Gaussian process surrogate model is then fitted to this training data in order to interpolate the model in parameter space for inference. The number and selection of training points significantly affects the performance of the surrogate model.

Interestingly, there is a strong connection between the unscented transform and Gaussian process surrogate models. It can be shown that sigma points for the UT are optimal training points for a Gaussian process surrogate model in the sense of minimizing the variance of the expected value of the posterior distribution (Särkkä et al., 2016). While the details are technical, one can think of UT sigma points as optimal training points given that (i) the prior distribution is a multivariate Gaussian and (ii) the Gaussian process surrogate model has a polynomial covariance function.

5.1 Modeling conclusions

Our work follows a number of previous observationally constrained paleo-ice sheet modeling studies (e.g., Tarasov and Peltier, 2002; Lecavalier et al., 2014; Calov et al., 2015). Perhaps most relevant to this work is Lecavalier et al. (2014), who model the deglaciation of Greenland from the Last Glacial Maximum (LGM) using a 3D thermomechanically coupled ice sheet model. Model runs in Lecavalier et al. (2014) are informed by constraints on relative sea level, ice core thinning, and LGM ice sheet extent.

In contrast to previous modeling studies, the computational efficiency of the flow line model outlined in Brinkerhoff et al. (2017) makes time-dependent data-assimilation, sensitivity testing, and robust uncertainty estimation tractable. Sensitivity testing indicates that estimated precipitation is insensitive to parameter uncertainties in the PDD and ice dynamics models. This conclusion supports earlier findings showing that modeled Holocene retreat in land-terminating sectors of the GrIS is more sensitive to surface mass balance than other factors like the flow law or basal sliding (Cuzzone et al., 2019; Lecavalier et al., 2014).

A drawback of our modeling approach is that we cannot account for inherent map plane effects such as changes in ice flow direction, or convergent or divergent flow in or out of the assumed flow path. These factors likely contribute to small discrepancies in estimated precipitation anomalies between the northern and southern flow lines (Fig. 6). Beyond computational efficiency, there are number of reasons why a flow line model is appropriate for our purposes. The modern flow field near Kangerlussuaq can be characterized by relatively simple east-to-west flow, and there is an absence of strongly convergent flow into outlet glaciers. Given the low bedrock relief in the region, and the surface-mass-balance-dominated retreat pattern (Van Tatenhove et al., 1996; Cuzzone et al., 2019), it is reasonable to assume that the flow regime was similar during the Holocene. In this work, we do not treat temperature as a random variable with its own covariance structure. However, differences between the Buizert et al. (2018) and Dahl-Jensen et al. (1998) inversions indicate that temperature is the dominant source of uncertainty in ΔP. This result underscores the importance of generating improved, regionally specific, temperature reconstructions constrained by proxy records. More regionally specific estimates of temperature would help to decrease uncertainty in the estimated precipitation history.

Despite lingering uncertainties, our modeling results indicate that the Holocene thermal maximum was accompanied by elevated snowfall, which slowed ice retreat in the Kangerlussuaq region of the GrIS. Inversions conducted using both the Buizert et al. (2018) and Dahl-Jensen et al. (1998) temperature reconstructions show average HTM snowfall around 35 % higher than modern. More specifically, in the Dahl-Jensen et al. (1998) inversion, there is a clear trend between HTM warming and increased accumulation, which could be interpreted as transient increase in precipitation due to reduced Arctic sea ice cover or changes in atmospheric circulation.

Appendix A: A higher-order method for estimating covariance

Estimating the mean and covariance of the posterior distribution requires approximating Gaussian-weighted expectation integrals of the form

$$\begin{array}{}\text{(A1)}& E\left[\mathcal{F}\right(\mathit{x}\left)\right]=\underset{{\mathbb{R}}^{\mathrm{n}}}{\int }\mathcal{F}\left(\mathit{x}\right)\mathcal{N}\left(\mathit{x}\right|\mathbf{0},I)\text{d}\mathit{x},\end{array}$$

via numerical integration rules, also called cubature rules, which approximate the expectation integral as a weighted sum

$$\begin{array}{}\text{(A2)}& E\left[\mathcal{F}\right(\mathit{x}\left)\right]\approx Q\left[\mathcal{F}\right(\mathit{x}\left)\right]=\sum _{i=\mathrm{1}}^N{w}_i\mathcal{F}\left({\mathit{\chi }}_i\right).\end{array}$$

Here, the notation $\mathcal{N}\left(\mathit{x}\right|\mathbf{0},I$) refers to a Gaussian probability density function evaluated at the point x. General Gaussian weight functions $\mathcal{N}\left(\mathit{x}\right|{\mathit{x}}_{\mathrm{0}},{\mathbf{P}}_x$) are handled by changing variables. Letting $\sqrt{{\mathbf{P}}_x}$ be a matrix square root of the covariance matrix, we have

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& \underset{{\mathbb{R}}^{\mathrm{n}}}{\int }\mathcal{F}\left(\mathit{x}\right)\mathcal{N}\left(\mathit{x}\right|{\mathit{x}}_{\mathrm{0}},{\mathbf{P}}_x)\mathrm{d}\mathit{x}\\ & =\underset{{\mathbb{R}}^{\mathrm{n}}}{\int }\mathcal{F}({\mathit{x}}_{\mathrm{0}}+\sqrt{{\mathbf{P}}_x}\mathit{\xi })\mathcal{N}\left(\mathit{\xi }\right|\mathbf{0},I)\mathrm{d}\mathit{\xi }\end{array}.\end{array}$$

Cubature rules, including the unscented transform, are constructed to exactly integrate polynomial functions ℱ(x) up to a certain degree d. Suppose that $\mathit{x}=[x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_n{]}^T$ is a point in ℝⁿ. A monomial of degree d refers to a function $x_{\mathrm{1}}^{i_{\mathrm{1}}}{x}_{\mathrm{2}}^{i_{\mathrm{2}}}\mathrm{\dots }{x}_n^{i_n}$, where the exponents are nonnegative integers that sum to d. A polynomial of degree d is a linear combination of monomials with highest degree d.

Li et al. (2017) describe a fifth-order cubature rule using fully symmetric sets of sigma points. A set $X=\mathit{\{}{\mathit{\chi }}_{\mathbf{1}},{\mathit{\chi }}_{\mathbf{2}},\mathrm{\dots },{\mathit{\chi }}_{\mathit{N}}\mathit{\}}$ is fully symmetric if it is closed under the operations of coordinate position and sign permutations. Their cubature rule has the form

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}Q\left[\mathcal{F}\right(\mathit{x}\left)\right]& =w_{\mathrm{1}}\mathcal{F}(\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{0})\\ & +w_{\mathrm{2}}\sum _{\text{full sym}}\mathcal{F}(\mathit{\lambda },\mathrm{0},\mathrm{\dots },\mathrm{0})\\ & +w_{\mathrm{3}}\sum _{\text{full sym}}\mathcal{F}(\mathit{\lambda },\mathit{\lambda },\mathrm{0},\mathrm{\dots },\mathrm{0}).\end{array}\end{array}$$

The notation ${\sum }_{\text{full sym}}\mathcal{F}(\cdot$) refers to a sum of the function ℱ evaluated at all points in the fully symmetric set generated by the given point.

Due to the symmetry of the sigma points and the Gaussian weight function, all moments (that is, integrals of Gaussian-weighted monomial functions) containing an odd order exponent are automatically satisfied. Exploiting this fact and the symmetries of the sigma points, it can be shown that satisfying the remaining moment constraint equations up the fifth order reduces to solving the following system of four equations in four unknowns $w_{\mathrm{1}},w_{\mathrm{2}},w_{\mathrm{3}}$ and λ

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& E\left[\mathrm{1}\right]=w_{\mathrm{1}}+\mathrm{2}n{w}_{\mathrm{2}}+\mathrm{2}n(n-\mathrm{1})w_{\mathrm{3}},\\ & E\left[x_i^{\mathrm{2}}\right]=\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{w}_{\mathrm{2}}+\mathrm{4}(n-\mathrm{1}){\mathit{\lambda }}^{\mathrm{2}}{w}_{\mathrm{3}},\\ & E\left[x_i^{\mathrm{4}}\right]=\mathrm{2}{\mathit{\lambda }}^{\mathrm{4}}{w}_{\mathrm{2}}+\mathrm{4}(n-\mathrm{1}){\mathit{\lambda }}^{\mathrm{4}}{w}_{\mathrm{3}},\\ & E\left[x_i^{\mathrm{2}}{x}_j^{\mathrm{2}}\right]=\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}{w}_{\mathrm{3}}.\end{array}\end{array}$$

By slightly modifying this cubature rule

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}Q\left[\mathcal{F}\right(\mathit{x}\left)\right]& =w_{\mathrm{1}}\mathcal{F}(\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{0}),\\ & +w_{\mathrm{2}}\sum _{\text{full sym}}\mathcal{F}({\mathit{\lambda }}_{\mathrm{1}},\mathrm{0},\mathrm{\dots },\mathrm{0}),\\ & +w_{\mathrm{3}}\sum _{\text{full sym}}\mathcal{F}({\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{2}},\mathrm{0},\mathrm{\dots },\mathrm{0}),\end{array}\end{array}$$

we introduce a new free parameter λ₂ that allows scaling of the sigma points about the mean. The new moment constraint equations become

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& E\left[\mathrm{1}\right]=w_{\mathrm{1}}+\mathrm{2}n{w}_{\mathrm{2}}+\mathrm{2}n(n-\mathrm{1})w_{\mathrm{3}},\\ & E\left[x_i^{\mathrm{2}}\right]=\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}{w}_{\mathrm{2}}+\mathrm{4}(n-\mathrm{1}){\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}{w}_{\mathrm{3}},\\ & E\left[x_i^{\mathrm{4}}\right]=\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{4}}{w}_{\mathrm{2}}+\mathrm{4}(n-\mathrm{1}){\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{w}_{\mathrm{3}},\\ & E\left[x_i^{\mathrm{2}}{x}_j^{\mathrm{2}}\right]=\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{w}_{\mathrm{3}}.\end{array}\end{array}$$

Using E[1]=1, $E\left[x_i^{\mathrm{2}}\right]=\mathrm{1}$, $E\left[x_i^{\mathrm{4}}\right]=\mathrm{3}$, and $E\left[x_i^{\mathrm{2}}{x}_j^{\mathrm{2}}\right]=\mathrm{1}$ we obtain

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& {\mathit{\lambda }}_{\mathrm{1}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}\sqrt{n-\mathrm{4}}}{n-{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{1}}},\\ & w_{\mathrm{2}}={\displaystyle \frac{\mathrm{4}-n}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{4}}}},\\ & w_{\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}},\\ & w_{\mathrm{1}}=\mathrm{1}-\mathrm{2}n{w}_{\mathrm{2}}-\mathrm{2}n(n-\mathrm{1})w_{\mathrm{3}},\end{array}\end{array}$$

with n>4 and $n-{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{1}\ne \mathrm{0}$. A drawback of the original cubature rule, as well our modified version here, is that it requires negative weights, which can lead to numerical instability.