1. Introduction
The majority of the Antarctic ice sheet is buttressed by floating ice shelves (Bindschadler and others, Reference Bindschadler2011) which provide large contact areas for ice–ocean interactions. Approximately half of the ice shelves’ total mass loss is attributed to ocean-induced melting at the underside of ice shelves (Depoorter and others, Reference Depoorter2013), and its spatiotemporal variability imprints ice flow dynamics farther upstream (Reese and others, Reference Reese, Gudmundsson, Levermann and Winkelmann2017; Gudmundsson and others, Reference Gudmundsson, Paolo, Adusumilli and Fricker2019). Consequently, ice flow and ocean models need to be coupled for future projections; frameworks (Goldberg and others, Reference Goldberg, Gourmelen, Kimura, Millan and Snow2019; Gladstone and others, Reference Gladstone2021), parameterizations (Burgard and others, Reference Burgard, Jourdain, Reese, Jenkins and Mathiot2022; Goldberg and Holland, Reference Goldberg and Holland2022) and benchmarks (Asay-Davis and others, Reference Asay-Davis2016) for this task have been developed. Similarly, the local snow accumulation is influenced by atmospheric conditions and is crucial in determining ice shelf thickness (Winkelmann and others, Reference Winkelmann, Levermann, Martin and Frieler2012). As a result, ice flow models are also coupled to climate models for future projections (Goelzer and others, Reference Goelzer, Huybrechts, Loutre and Fichefet2016; Pattyn and others, Reference Pattyn, Favier, Sun and Durand2017). It is crucial to confront ice flow models with observations to validate them and investigate their ability to explain observed phenomena. Here, we present a new method that infers surface accumulation (also known as ‘surface mass balance’ (Lenaerts and others, Reference Lenaerts, Medley, van den Broeke and Wouters2019)) and basal melt rates (collectively, the mass-balance parameters) from the ice shelves’ internal stratigraphy, which can be routinely mapped by radio-echo sounding.
Typically, surface accumulation is the more accessible mass-balance parameter (Eisen and others, Reference Eisen2008); it can be measured in situ using stake farms and can also be derived from multiple firn cores (Lenaerts and others, Reference Lenaerts, Medley, van den Broeke and Wouters2019). Many of these observations validate atmospheric models such as RACMO (van Wessem and others, Reference van Wessem2018) and MAR (Gallée and Schayes, Reference Gallée and Schayes1994; Agosta and others, Reference Agosta2019), which estimate surface accumulation on 
$35\,\mathrm{km}$ grids (Lenaerts and others, Reference Lenaerts, Medley, van den Broeke and Wouters2019) (with few locations being estimated at a higher resolution of 
$5.5\,\mathrm{km}$). Estimating the basal melt is more challenging and is typically dependent on knowledge of surface accumulation. For example, estimates of surface accumulation have been used along with mass conservation arguments to estimate basal melt (Neckel and others, Reference Neckel, Drews, Rack and Steinhage2012; Depoorter and others, Reference Depoorter2013; Berger and others, Reference Berger, Drews, Helm, Sun and Pattyn2017; Adusumilli and others, Reference Adusumilli, Fricker, Medley, Padman and Siegfried2020). These approaches have provided Antarctic-wide time series of the last few decades of basal melt rates (Adusumilli and others, Reference Adusumilli, Fricker, Medley, Padman and Siegfried2020). The spatial resolution is currently limited to the kilometer scale, which may miss fine grained processes occurring within ice shelf channels (Drews, Reference Drews2015; Marsh and others, Reference Marsh2016) or near basal terraces(Dutrieux and others, Reference Dutrieux2014). Measurements of basal melt which are independent of the surface accumulation are also available, but typically only on short temporal scales, for example, with time-lapse radar measurements of ice thickness change (Zeising and others, Reference Zeising, Steinhage, Nicholls, Corr, Stewart and Humbert2022). Using phase-coherent data acquisition, these measurements can disentangle the observed thickness change into strain thinning and basal melt (Nicholls and others, Reference Nicholls, Corr, Stewart, Lok, Brennan and Vaughan2015). This has provided much insights, e.g., in terms of relevant tidal (Sun and others, Reference Sun, Hattermann, Pattyn, Nicholls, Drews and Berger2019) and seasonal timescales (Vankova and Nicholls, Reference Vankova and Nicholls2022).
Here, we investigate to what extent the radar-imaged isochronal ice stratigraphy (Eisen and others, Reference Eisen, Nixdorf, Wilhelms and Miller2004) can provide additional information for inferring mass-balance parameters. On grounded ice, radar-imaged internal reflection horizons (IRHs) have been used in multiple ways, for example, to infer the surface accumulation history (Waddington and others, Reference Waddington, Neumann, Koutnik, Marshall and Morse2007; MacGregor and others, Reference MacGregor, Matsuoka, Koutnik, Waddington, Studinger and Winebrenner2009; Catania and others, Reference Catania, Hulbe and Conway2010; Steen-Larsen and others, Reference Steen-Larsen, Waddington and Koutnik2010; Wolovick and others, Reference Wolovick, Moore and Zhao2021; Theofilopoulos and Born, Reference Theofilopoulos and Born2023), velocity patterns of the ice flow (Eisen, Reference Eisen2008; Holschuh and others, Reference Holschuh, Parizek, Alley and Anandakrishnan2017), ice-rise evolution (Drews and others, Reference Drews, Matsuoka, Martín, Callens, Bergeot and Pattyn2015; Henry and others, Reference Henry2023) or large-scale model calibration (Leysinger Vieli and others, Reference Leysinger Vieli, Hindmarsh, Siegert and Bo2011; Sutter and others, Reference Sutter, Fischer and Eisen2021). On ice shelves, surface accumulation rate can also be derived from the radar-measured shallow stratigraphy (Pratap and others, Reference Pratap2022), but not from intermediate depths and below where the stratigraphy is also influenced by basal melt and ice flow. The stratigraphy of ice shelves differs for various combinations of surface accumulation and basal melt rates (Visnjevic and others, Reference Visnjevic2022). This suggests that given an ice flow model of the internal stratigraphy that accounts for the surface accumulation and basal melt rates, we can use observed IRHs to recover the surface accumulation and basal melt rate histories (Fig. 1). Thus, our goal is to solve the inverse problem of inferring the surface accumulation and basal melt rates that can explain the observed IRHs under the physical constraints of the ice flow model.
Figure 1. Estimation of mass-balance parameters from a steady-state ice shelf with two methods. The Eulerian Mass Budget method (left) detects the difference of surface accumulation and basal melt within two flux gates (blue vertical lines) by considering flux divergence 
$\nabla\cdot(\bf{vh})$. Often, the basal melt rates 
$\dot b$ are inferred assuming that the surface accumulation (
$\dot{a}_{\text{obs}}$) is known. In the internal reflection horizon (IRH) method, we are given information on the internal stratigraphy of the shelf. This information is used to separate the known total mass balance into individual estimates of surface accumulation and basal melt (
$\dot{a}_{\text{avg}},\dot{b}_{\text{avg}}$ respectively). These estimates correspond to the time-averaged value over the age of the IRH to the present. The inset plots show different surface accumulation and basal melt parameterizations which give rise to the same total mass balance and overall shape of the ice shelf, but different internal stratigraphy.
Inverse problems, also known as inversion, data assimilation or inference problems in the literature, denote the task of finding the model parameters that are compatible both with empirical observations and prior knowledge. This problem is widespread in the geosciences, e.g., in hydrogeology (Linde and others, Reference Linde, Renard, Mukerji and Caers2015), seismology (Symes, Reference Symes2009) or in climate science (Tebaldi and Sansó, Reference Tebaldi and Sansó2008). Bayesian inference provides a powerful framework for solving inference problems, but conventional Bayesian approaches are restricted to models for which the so-called ‘likelihood function’ is tractable. A tractable likelihood function is one that can be efficiently evaluated (see Appendix C.1 for examples of tractable and intractable likelihood functions). However, this is not the case in our setup. We therefore use simulation-based inference (SBI, Papamakarios and Murray (Reference Papamakarios and Murray2016); Lueckmann and others (Reference Lueckmann, Goncalves, Bassetto, Öcal, Nonnenmacher and Macke2017); Cranmer and others (Reference Cranmer, Brehmer and Louppe2020)) to solve the inverse problem presented in this work. In SBI, we evaluate the forward model under different values of the model’s parameters from a prior distribution. We use the resulting simulated dataset to train a neural network that performs conditional density estimation. In the neural posterior estimation (NPE) variant, the network approximates the Bayesian posterior distribution. A key advantage of NPE is the amortization of simulation cost. An amortized inference framework is one that, once trained, can be instantly applied to find the posterior distribution for any new observation without requiring more simulations or training. Importantly, SBI does not require the forward model to be differentiable and can also work with ‘blackbox’ models. Therefore, our approach can be extended to a variety of preexisting forward models. To the authors’ knowledge, this work is the first application of SBI in glaciology, but we note that it has already been applied in other geoscientific disciplines such as geothermics (Omagbon and others, Reference Omagbon2021), hydrogeology (Allgeier and Cirpka, Reference Allgeier and Cirpka2023), hydrology (Hull and others, Reference Hull2022) and molecular ecology (Overcast and others, Reference Overcast2021).
In this study, we consider steady-state ice shelves and IRHs in the local meteoric ice (LMI) body of ice shelves (Das and others, Reference Das2020). This work is a test case for inferring atmospheric and oceanographic boundary conditions from the ice stratigraphy with a novel inference technique that provides uncertainty estimates. Our approach can be transferred to other ice flow regimes (e.g. flank flow on grounded ice) where similar scientific questions can be explored. Our approach can similarly be adapted to ice shelves exhibiting marine ice formation. Moreover, the isochronal stratigraphy of ice shelves and ice sheets (including the neighboring ice rises) is currently the only archive of surface accumulation and basal melt over the past hundreds of years. Our approach is capable of testing this archive. Thus, this study provides one link between observational initiatives (such as AntArchitecture, Bingham and others, Reference Bingham2024) for Antarctica-wide internal stratigraphy datasets and the modeling community.
The paper is structured as follows: In Section 2, we describe our forward model of the internal stratigraphy of an ice shelf and introduce our inference approach. In Section 3, we detail the synthetic ice shelf construction. We also present the results of inferring the mass-balance parameters from this synthetic stratigraphy and compare the posterior distribution to a known ground truth. In Section 4, we describe the setting of the Ekström Ice Shelf (EIS) and the dataset of observed IRHs along the central flow line transect. We then provide the results of our inference framework and compare them to independent measurements of surface accumulation uniquely available in this location for the periods 1996–2005 and 2014–23. In Section 5, we interpret our results and evaluate our approach. We finally conclude and discuss future perspective in Section 6.
2. Methodology
2.1. Forward model
We denote spatially varying parameters as functions, e.g. 
$\dot{a}(x)$ or at times 
$\dot{a}$ for brevity, while bold-faced characters denote the discretized values of this function on a specified grid, e.g. 
$\dot{\mathbf{a}} = [\dot{a}(x_{1}),\ldots, \dot{a}(x_{n})]^{\top}$.
2.1.1. Ice flow model
We model ice shelves using the shallow shelf approximation (SSA) (Morland, Reference Morland1984). Throughout this study, we consider ice shelves in steady state. Consequently, the ice surface s, base f, thickness 
$h = s-f$ and velocity v are all fixed throughout our simulations. We assume plug flow for the ice shelf regime, meaning that the horizontal velocity profile does not change in the vertical direction z. These assumptions results in the mass-balance condition
\begin{equation}
\nabla \cdot (hv) = \dot{m} ,
\end{equation} where hv is the total mass flux, 
$\nabla \cdot$ is the divergence operator and 
$\dot{m} = \dot{a} -\dot{b}$ is the total mass-balance rate. Here we use the convention that the surface accumulation rate 
$\dot{a}$ is positive for mass gain of the ice shelf and the basal melt rate 
$\dot{b}$ is positive for mass loss. In this exploratory study, we focus on flow lines. We parameterize our domain such that x denotes the distance along the flow line, and vx now denotes the velocity parallel to the flow line. The two-dimensional (2-D) geometry is only valid for observations located on flow lines and in the absence of lateral compression and extension. While the former is approximately true in our case, the latter is unrealistic for most Antarctic ice shelves. To account for ice flux into or out of our modeling domain, we, therefore, include the ice flux component normal to the flow line as an additional, spatially variable term to the total mass-balance rate 
$\dot{m}$ (Appendix A). We test the validity of this approach in a 2-D synthetic example (Section 3) that includes a spatially variable total mass balance and lateral compression. For the real-world scenario, we estimate the the normal ice flux component from satellite velocities.
We seek to predict the steady-state internal stratigraphy for a given flow line and possible surface accumulation and basal melting rate profiles. We define the internal stratigraphy to be a set of isochronal layer elevations 
$\{e_{1}(x),\ldots e_{L}(x)\}$, with 
$f(x) \leq e_{1}(x)\leq e_{2}(x)\leq \cdots \leq e_{L}(x) \leq s(x)$. One approach to calculate the internal stratigraphy uses the SSA expression for the vertical component of the velocity (Greve and Blatter, Reference Greve and Blatter2009) to have a fully specified velocity field. This can then be used to calculate the age field 
$\mathcal{A}(x,z)$ of the shelf. Contours of constant age (isochrones) then define the internal stratigraphy. However, these methods suffer from numerical diffusion and can be computationally expensive (Visnjevic and others, Reference Visnjevic2022).
The computational efficiency of the forward model is crucial for our inference method, as we need to evaluate the forward model many times. As a result, we opt instead to use an implementation of the tracer method (Born, Reference Born2017; Born and Robinson, Reference Born and Robinson2021). The model is seeded with vertical segments each with a thickness profile 
$\{h_{1}(x),\ldots, h_{L}(x)\}$, such that the sum matches the ice geometry 
$\sum_{l=1}^{L}h_{l}(x) = h(x)$. The horizontal velocity 
$v_{x}(x)$ is used to advect mass within segments and to thin or thicken the segments as a function of the prescribed strain rates. The accumulation 
$\dot{a}(x)$ and melt 
$\dot{b}(x)$ rates are used to add new segments or take away mass from the two boundary segments at the top and bottom of the shelf respectively. The (isochronal) layer elevations are then the boundaries between our modeled segments. We use the convention that el corresponds to the top of segment l, which can be calculated using the cumulative thicknesses of the segments below,
\begin{equation}
e_{l}(x) = f(x) + \sum_{l'=1}^{l} h_{l'}(x).
\end{equation}In our simulations, we used a high temporal resolution of one isochronal layer per year. Despite the high resolution, the layer tracing method allows for determining the internal stratigraphy in a computationally efficient manner. For the domains and timescales considered in our study, the complete forward model can be evaluated on the order of 60 s on a single CPU core, enabling the application of SBI methods (see Appendix E for details).
To uniquely determine the layer thicknesses in such a scheme, we need to specify the boundary conditions on the layer thicknesses hl at the inflow boundary x = 0 (here corresponding to the grounding line). The true boundary conditions are typically not known. However, the stratigraphy in a large part of the domain is still independent of the boundary conditions. This zone corresponds to the LMI body of ice shelves (Das and others, Reference Das2020). When inferring from observed stratigraphy data, we use only data within the LMI body. We detail our model of the LMI body in Appendix A.
2.1.2. Noise model
The ice flow model predicts isochronal layers with varying depth over spatial scales of kilometers. Observed IRHs, however, also show variability on sub-kilometer scales. This systematic model-data misfit is caused by errors in input datasets (such as surface velocity, geometry), coarse resolution of the forward model and omission of higher order processes that are not included in the forward model, such as the effect of rheology. For inference, it is important that the predicted isochrones have consistent statistical properties with the observed IRHs. This is achieved by the definition of an appropriate noise model.
The ice flow model predicts isochronal layer elevations 
$\{\mathbf{e}_{1}(\mathbf{x}),\ldots,\mathbf{e}_{L}(\mathbf{x})\}$ on a fixed grid 
$\mathbf{x}\in {\mathbb{R}}^{N}$ where N is the number of grid points. Guided by empirical observations, the noise model should have the property that the errors of different between modeled layers l at different depths are spatially correlated and amplified for deeper layers. We, therefore, define a layer-wise noise model as the product of an x-dependent baseline noise function and a z-dependent vertical amplification factor. More precisely, the additive noise 
$\boldsymbol{\delta}_{l}\in\mathbb{R}^{N}$ of layer l is defined as
\begin{equation}
\boldsymbol{\delta}_{l} = \boldsymbol{\epsilon} \odot \mathbf{T}(\mathbf{e}_{l}),
\end{equation} where 
$\boldsymbol{\epsilon} = [\epsilon_{1},...,\epsilon_{N}]^{\top}$ is a x-dependent noise profile, which is shared for all layers, 
$\mathbf{T}(\mathbf{e}_{l}) = [T(e_{l,1}),\ldots,T(e_{l,N})]^{\top}$ is a deterministic function of elevation (increasing with depth), and 
$\odot$ denotes an element-wise product. The vertical scaling 
$\mathbf{T}(\cdot)$ mimics uncertainties in the travel time-to-depth conversion which depend on the density 
$\rho(z)$. Here, this is done using 
$\rho(z)$ as in Drews and others Reference Drews(2016) and an empirical density-permittivity relation (Looyenga, Reference Looyenga1965) to calculate the radio-wave speed c(z). This results in the factor
\begin{equation}
T(z) = \int_{z}^{s}\frac{dz'}{c(z')},
\end{equation}which we then discretize on the set of layer elevations.
The sub-kilometer variability of the observed IRHs are modeled with power spectral densities ϵ:
\begin{equation}
\boldsymbol{\epsilon} = A_{\boldsymbol{\epsilon}}\sum_{n=1}^{N} \sqrt{\exp^{\beta_{n}}}\cos(2\pi\omega_{n} + \chi_{n}),
\end{equation} where the log power spectral densities βn and offsets χn are randomly sampled from normal and uniform distributions respectively: 
$\beta_{n}\sim\mathcal{N}(\mu_{\beta_{n}},\sigma_{\beta_{n}}^{2})$ and 
$\chi_{n}\sim U([-\pi,\pi))$. The frequencies ωn are the corresponding Fourier frequencies of the simulation grid x and 
$A_{\boldsymbol{\epsilon}}$ is a global scale factor (set to 
$4\times 10^{-10}$). In the synthetic ice shelf (Section 3), we define the distribution of the log power densities using 
$\sigma_{\beta_{n}}^{2} = 0.5$ and
\begin{equation}
\mu_{\beta_{n}} = -8\left( 1-\exp^{-200\omega_{n}}\right).
\end{equation} For EIS, the distribution means 
$\mu_{\beta_{n}}$ and variances 
$\sigma_{\beta_{n}}^{2}$ were calibrated given the observed IRHs on a separate set of calibration simulations (full details in Appendix B). We emphasize that this representation of the noise model is a choice—we define a mathematical model of the mismatch, rather than model a physical effect directly. Thus, other choices are possible. We choose this representation of the noise model for its flexibility and interpretability.
By combining the ice flow model with the empirically guided noise model, we have arrived at a physically motivated forward model to sample a plausible observed internal stratigraphy of an ice shelf from the mass-balance rate parameters 
$\dot{a}$ and 
$\dot{b}$.
2.2. Inference
Having established the forward model, we arrive at the inverse problem of finding the surface accumulation and basal melt rates that best explain the observed internal stratigraphy. We use Bayes theorem with model parameters θ and outcomes X:
\begin{equation}
p(\theta|X) = \frac{p(X|\theta)p(\theta)}{p(X)} .
\end{equation} Here, 
$p(\theta|X)$ is the posterior distribution of the parameters given a particular outcome X, 
$p(X|\theta)$ is the likelihood function of the model, 
$p(\theta)$ is the prior distribution encoding our existing knowledge on the plausible values of θ and p(X) is the model evidence. The goal of Bayesian inference is to find the posterior 
$p(\theta|X_{o})$, where Xo is observed data which has the same form as X, but is measured, instead of simulated.
2.2.1. Simulation-based inference
It is generally not possible to analytically solve for the Bayesian posterior distribution (Eqn (7)), as the evidence term p(X) cannot be computed. Approximate methods exist to solve Eqn (7) using knowledge of only the likelihood function and prior distribution. In this work, we deploy SBI, an approximate Bayesian inference and likelihood-free approach, using only samples from our forward model. In SBI, we use artificial neural networks to approximate conditional probability distributions. While there exist different variants of SBI which target either the likelihood 
$p(X|\theta)$ or the likelihood ratio (see Cranmer and others Reference Cranmer, Brehmer and Louppe(2020) for an overview), we focus on NPE, which approximates the posterior distribution directly (Papamakarios and Murray, Reference Papamakarios and Murray2016; Greenberg and others, Reference Greenberg, Nonnenmacher and Macke2019).
In NPE, we generate a training dataset 
$\{(\theta_{k},X_{k})\}_{k=1}^{K}$ (Fig. 2) by sampling parameters from the prior 
$\theta_{k}\sim p(\theta)$ and sampling from the forward model 
$X_{k} \sim p(X|\theta_{k})$. To approximate the posterior distribution, a variational family of distributions 
$q_{\phi}(\theta|X)$ is typically defined in terms of a neural network with learnable weights ϕ. We represent qϕ as a normalizing flow (Durkan and others, Reference Durkan, Bekasov, Murray and Papamakarios2019; Kobyzev and others, Reference Kobyzev, Prince and Brubaker2019; Papamakarios and others, Reference Papamakarios, Nalisnick, Rezende, Mohamed and Lakshminarayanan2019a). Normalizing flows are flexible generative models, which, once trained, can be used either to sample or evaluate the (log-) probability density function of the conditional distribution 
$q_\phi(\theta|X)$, for any outcome X in the support of the training dataset. We provide a brief description of normalizing flows in Appendix C.2 and refer the reader to Papamakarios and others Reference Papamakarios, Nalisnick, Rezende, Mohamed and Lakshminarayanan(2019a) for a review.
Figure 2. Simulation-based inference workflow. SBI has two primary phases: training and evaluation. In the training phase, accumulation rates are randomly sampled from a prior distribution, the corresponding basal melt rates are obtained using total mass balance, and the resulting internal stratigraphy is calculated using the forward model. These simulations from the prior are used to train a neural network which parameterizes conditional distributions. In the evaluation phase, the trained network is conditioned on the observed IRH and outputs the Bayesian posterior distribution over the parameters (without any additional calls to the forward model).
In NPE, the neural network is trained by minimizing the expected negative log-probability
\begin{equation}
\mathcal{L}(\phi) = \mathbb{E}_{\theta_k\sim p(\theta),X_k\sim p(X|\theta_k)}[ -\log q_{\phi}(\theta_k|X_k)]
\end{equation}on the training dataset. More intuitively, this loss seeks to maximize the probability assigned to the training data. It can be trivially shown that minimizing this loss is equivalent to minimizing the (forward) Kullback–Leibler (KL) divergence between the variational distribution and the true posterior distribution (see C.3).
It has been shown that, if there exists a set of weights ϕ such that 
$q_\phi(\theta|X)$ is the true posterior distribution, and in the limit of infinite training samples 
$K\to\infty$, the minimum of the loss in Eqn (8) is reached when 
$q_{\phi}(\theta|X) = p(\theta|X)$ for all X—i.e. when our estimated distribution matches the true posterior (see Proposition 1 of Papamakarios and Murray Reference Papamakarios and Murray2016 for full statement and proof).
We additionally make use of an embedding network, which are commonly used in SBI workflows to improve performance. Embedding networks learn summary statistics Y(X), which are lower-dimensional representations of the outcomes X. Using the embedding Y(X) as an input to the normalizing flow instead of X itself reduces the model complexity. The embedding network is trained jointly with the normalizing flow. In our setting, Xk are spatially varying IRH elevations, and so we choose a 1-dimensional (1-D) convolutional neural network as our embedding net, resulting in 50-dimensional embeddings on which the posterior network is conditioned (full details in Appendix D). Throughout this work, we use the sbi package for Python (Tejero-Cantero and others, Reference Tejero-Cantero2020) to perform inference.
2.2.2. Definitions of model parameters, outputs and observations
We define 
$\theta = \dot{\mathbf{a}}_{\text{inf}}= [\dot{a}_{i_{1}},\ldots,\dot{a}_{i_{J}}]^{\top}$, the values of the surface accumulation rate on a discretized grid 
$\tilde{\mathbf{x}}$. In our experiment, we choose the number of inference grid points J = 50 as a compromise between computational complexity of the inference problem while still inferring accumulation rate at a high resolution of ∼2.5 km. This is smaller than the discretized grid x we use for our simulations, which has 500 gridpoints in our experiments. In practice, we take 
$\tilde{\mathbf{x}}$ to be a regularly spaced subset of x, so that 
$\dot{\mathbf{a}}_{\text{inf}}$ can also be taken as a subset of a. However, 
$\tilde{\mathbf{x}}$ can be any discretization of the flow line and need not be a subset of x. Furthermore, despite defining θ to only represent the surface accumulation, any inference of the surface accumulation rate automatically extends to inference of the basal melt rates. This is because for any probability distribution 
$q(\dot{\mathbf{a}}_{\text{inf}})$, the total mass-balance relationship implies that 
$\dot{\mathbf{b}}_{\text{inf}} \sim q(\dot{\mathbf{a}}_{\text{inf}} - \dot{\mathbf{m}}_{\text{inf}})$, where 
$\dot{\mathbf{b}}_{\text{inf}},\dot{\mathbf{m}}_{\text{inf}}$ are the respective discretizations of 
$\dot{b},\dot{m}$ onto 
$\tilde{\mathbf{x}}$.
We now turn to describing the observation, Xo, and forward model outcomes Xk. The observed data are a set of different IRHs, 
$\{\mathbf{e}_{m}(\mathbf{x})\}_{m=1}^{M}$, where 
$e_{m}(x_{i})$ is the elevation of the 
$m{\text{th}}$ IRH in our dataset at grid position i. The IRH elevations need not and typically are not observed at the same locations as the simulation gridpoints; and so we first interpolate the IRH elevations onto the simulation grid x using linear spline interpolation (as implemented in Scipy (Virtanen and others, Reference Virtanen2020)). Therefore, we assume 
$\{\mathbf{e}_{m}(\mathbf{x})\}_{m=1}^{M}$ is already defined on x. One reasonable choice is to define the observation Xo as the entire set of all measured IRH elevations. However, in our work, we choose to separately infer the mass balance from each IRH in our observed dataset. This choice has two advantages: first, ordering IRHs by depth also corresponds to their reverse age order, with the oldest IRHs being the deepest. Thus, inferring the surface accumulation and basal melt rates for deeper IRHs corresponds to inferring the average rates over longer periods of time. By comparing the inferred mass-balance parameters obtained with different IRHs, we can reason whether or not our steady-state assumption is valid. The second advantage is practical—we seek a consistent representation of the observations that can be applied across ice shelves. Given a different ice shelf, there will be a different number of IRHs at different depths. Therefore, the embedding net for these data will have to have a different architecture for each ice shelf. In our representation, the embedding net can always be a 1-D convolutional net, as the observations are always 1-D vectors.
Thus, given a dataset of M observed IRHs, we have M inference problems to solve, where each observation corresponds to one IRH. It is, therefore, reasonable to take one isochronal layer of the simulated stratigraphy as the output of the forward model. For the 
$m{\text{th}}$ inference problem, we define the outcome of the forward model as the isochronal layer el that is closest to IRH m (in the mean square sense). More precisely, for inference problem m and simulation k, we define the observation of the forward model to be 
$X_{k}^{m} = \mathbf{e}_{l^{*}}(\mathbf{x}_{i\geqslant i(m)})$, where
\begin{equation}
l^{*} = \operatorname{\rm{arg\,min}}\limits_{l} ||\mathbf{e}_{l}(\mathbf{x}_{i\geqslant i(m)})-\mathbf{e}_{m}(\mathbf{x}_{i\geqslant i(m)})||_{2}^{2}.
\end{equation} Here, i(m) is the index of the boundary of the LMI body for IRH 
$\mathbf{e}_{m}(\mathbf{x})$. For 
$i \lt i(m)$, the IRH 
$e_m(x_i)$ is outside the LMI body and for 
$i\geqslant i(m)$ within the LMI body (see Appendix B for details). We further define 
$\mathbf{x}_{i\geqslant i(m)} = [x_{i(m)}, x_{i(m)+1}, ..., x_{N}]^{\top}$ as the restriction of the gridpoints x to within the LMI body of 
$\mathbf{e}_{m}(\mathbf{x})$. We correspondingly set the observation for IRH m to 
$X_{o}^{m} = \mathbf{e}_{m}(\mathbf{x}_{i\geqslant i(m)})$. Our choice to select the simulated isochronal layer that most closely matches the IRH is due to the true age of the IRH being unknown. This introduces degeneracy into the forward model—two simulations with different surface accumulation and basal melt rate parameterizations can produce isochronal layers with a similar geometry but different ages. It is, therefore, important to define the prior distribution appropriately, which we do in the following section.
2.2.3. Choice of prior distribution
We aim to approximate the posterior distribution
\begin{equation}
p(\dot{\mathbf{a}}| X_{o}^{m}) \propto p(X_{o}^{m}|\dot{\mathbf{a}})p(\dot{\mathbf{a}}).
\end{equation} The likelihood 
$p(X_{o}^{m}|\dot{\mathbf{a}})$ is not tractable but can be sampled from using the forward model. To specify the prior, we use the long-term snow accumulation observations of the Neumayer stations (Wesche and others, Reference Wesche2016; Wesche and Regnery, Reference Wesche and Regnery2022) over more than 30 years to define an empirically motivated prior for EIS, which we also use for the synthetic ice shelf. We first assume that localized surface melt (
$\dot{a} \lt 0$) is possible, but rare. We also observe that average rate of accumulation is ∼0.5 
$\mathrm{m}\,\mathrm{a}^{-1}$, and that the accumulation rate is almost everywhere under 
$2\,\mathrm{m}\,\mathrm{a}^{-1}$. Finally, we take the accumulation rate to vary smoothly in space. We define a prior distribution that satisfies these criteria, while still allowing for a broad range of surface accumulation rate profiles. We define the following generative process for 
$\dot{\mathbf{a}}$: first, we draw a sample 
$\boldsymbol{\alpha} = [\alpha_{1},\ldots,\alpha_{N}]^{\top}$ from a Gaussian process with mean function µ = 0 and a Matérn kernel with a Matérn-ν of 2.5 and a length scale of 2500 m (Rasmussen and Williams, Reference Rasmussen and Williams2005). We then independently sample an offset 
$\mu_{\text{off}}\sim \mathcal{N}(0.5,0.25^{2})$ and scale 
$\sigma_{\text{sc}}\sim U(0.1,0.3)$ parameter. Finally, we set 
$\dot{\mathbf{a}} = \sigma_{\text{sc}}\dot{\boldsymbol{\alpha}} + \mu_{\text{off}}\mathbf{1}$. We inspect the implicit prior this defines over the basal melt rates in Section 4.3.
Defining the prior in this way is sufficiently expressive to capture numerous accumulation rate profiles, while also restricting the samples to conform to empirical knowledge. Additionally, the prior is shared for all M inference problems we have defined, and one evaluation of the forward model provides an observation 
$X_{k}^{m}$ for each of the inference problems. Thus, the same training dataset can be used for all posterior networks in our SBI approach, significantly reducing the computational costs.
3. Synthetic test case
Before we apply the presented workflow to EIS, we showcase its applicability in a synthetic test case in which all parameters are known.
3.1. Configuration of shelf and flow line
We test our workflow on a 2-D flow tube geometry from which we extract a flow line to infer the prescribed surface and basal accumulation rates as done later in the case of EIS. The flow tube is modeled using icepack (Shapero and others, Reference Shapero, Badgeley, Hoffman and Joughin2021) on a grid 
$L_{x}=125\,\mathrm{km}\times L_{y}=10\,\mathrm{km}$, with the along-flow direction x and across-flow direction y. We prescribe a Dirichlet boundary condition at the inflow and lateral boundaries, with a constant thickness of h 0, and a constant along-flow velocity of 
$v_{0x}$. The outflow boundary is set to be a static calving front. We initialize with a zero centered, longitudinally symmetric across-flow velocity 
$v_{0y}$ on the lateral boundaries, resulting in a flow field that has convergence (i.e. mass input) on the center flow line. We prescribe a spatially variable total mass balance 
$\dot{m}$: In our experiments, we set 
$v_{ox} = 100\,\mathrm{m}\,\mathrm{a}^{-1},\, v_{0y}=\pm 20\,\mathrm{m}\,\mathrm{a}^{-1}$ at y = 0 and 
$y=L_y$, respectively,
\begin{equation}
\dot{m} = -0.6 - 0.05\frac{x}{L_x} +0.3\exp\left(-\left(\frac{x-0.7L_x}{0.1L_x}\right)^{2}\right).
\end{equation} We let the geometry evolve under the SSA approximation until steady state is reached. From the steady-state ice shelf, we choose a discretization of the central flow line, x, and extract the relevant variables along this flow line to define the internal stratigraphy model (Fig. 3). The numerical values for additional parameters for the spin-up are given in Appendix D. The variables we need are the surface s and base f elevations, the along-flow velocities vx, and the along- and across-flow flux divergences 
$d(\mathbf{v}_{x}\mathbf{h})/dx$, 
$d(\mathbf{v}_{y}\mathbf{h})/dy$. These define the total mass balance, since:
\begin{equation}
\dot{\mathbf{a}}-\dot{\mathbf{b}} = \frac{d(\mathbf{v}_{x}\mathbf{h})}{dx} + \frac{d(\mathbf{v}_{y}\mathbf{h})}{dy}.
\end{equation}
Figure 3. Two-dimensional flow tube domain setup for the synthetic example. Map view of the simulated ice shelf’s surface. Flow lines (gray lines) converge to the central flow line (red). Color indicates ice thickness. The input variables for the internal stratigraphy model are evaluated on the central flow line.
We then solve the inverse problem which accounts in this case for mass gain through lateral compression. We choose a random sample from the prior distribution as the ground truth, 
$\dot{\mathbf{a}}_{\text{GT}}$, from which 
$\dot{\mathbf{b}}_{\text{GT}}$ follows accordingly. The forward model is then sampled to obtain a set of ground truth layer elevations, 
$\mathbf{e}_{o}(\mathbf{x})$. From these layer elevations, we choose to perform inference for four layers of ages 50, 100 and 150, and 300 years (labeled 1–4 in ascending order of age). These ages roughly correspond to the range of ages of the IRHs that we expect to observe on ice shelves.
3.2. Inference results
We evaluate the trained neural posterior network on the ground truth isochronal layer of age 50 years. The inferred posterior mean for the surface accumulation rate parameter is close to the ground truth accumulation rate (Fig. 4a,c) with the ground truth lying within the 95% confidence intervals of the posterior distribution.
Figure 4. Prior and posterior (predictive) for the synthetic dataset. (a and c) Prior and posterior over surface accumulation and basal melt rates respectively for layer 1 of the synthetic ice shelf, of age 50 years. Solid line is the distribution mean, the shaded region represents the 5th and 95th percentiles. The ground truth (GT) parameters used to generate the reference isochronoal layer are shown in red. (b) Cross section of the ice shelf. Prior and posterior predictive distributions for the layer closest matching the ground truth isochronal layer. The vertical dashed line represents the LMI boundary for this isochronal layer. The posterior predictive reconstructs the observed layer with higher accuracy and lower uncertainty. The posterior predictive distribution of the age of the isochronal layer is 
$60^{+9}_{-12}$ years (meaning a median of 60 years, and 16th and 84th percentiles of 48 and 69 years, respectively). The average root-mean-square error (RMSE) relative to the GT isochronal layer is 
$3.9\,\mathrm{m}$ for the posterior predictive distribution and 
$11.5\,\mathrm{m}$ for the prior predictive distribution.
Next, we evaluate the forward model on samples from the posterior (and prior) distribution to get the respective predictive distributions. The prior predictive distribution (Fig. 4b, green) is the distribution over the internal layers generated by simulating the forward model with mass-balance parameters drawn from the prior distribution. The posterior predictive distribution (Fig. 4b, blue) is defined similarly by simulating with mass-balance parameters from the posterior distribution. The posterior predictive matches the ground truth isochronal layer with high fidelity. We calculate the RMSE of the predictive simulations relative to the ground truth layer elevations for 1000 simulations using prior and posterior samples. The average RMSE for the posterior predictive distribution is 
$3.9\,\mathrm{m}$, compared to 
$11.5\,\mathrm{m}$ for the prior predictive distribution. Uncertainties in the layer elevations are much smaller than those of the prior predictive distribution. This is in contrast to the posterior uncertainty over the mass-balance rates, which is still considerable. This showcases the importance of our uncertainty-aware approach: there is more than one parameterization of accumulation and basal melt rates that can lead to similar isochronal layers.
The posterior uncertainty is also reflected in the inferred age of the isochronal layer. We infer an age of 
$60^{+9}_{-12}$ years for this layer (meaning a median of 60 years, and 16th and 84th percentiles of 48 and 69 years, respectively). This value closely matches the age of the ground truth isochronal layer, which was not used during inference. Thus, we have produced an estimate of the age of the layer without requiring time intensive measurements such as ice cores. We report the posterior distributions for deeper synthetic layers in Appendix F.
Finally, while we do not use the isochronal layer elevations outside the LMI boundary for inference, we can still infer the surface accumulation and basal melt rates at these locations. This is because the values of surface accumulation rate and basal melt rate still affect the downstream isochronal layer elevations, and so the observed elevations in the LMI body still contain information about the mass-balance rates upstream of the LMI boundary. Thus, we are still able to infer the mass-balance rates for 
$x \lt 15\,\mathrm{km}$.
4. Ekström Ice Shelf
EIS is a medium-sized ice shelf located between the Sörasen and Halvfarryggen Ice Rises in Dronning Maud Land, East Antarctica (Fig. 5c). EIS makes for an appropriate study site since the steady-state assumption likely holds (Drews and others, Reference Drews, Martín, Steinhage and Eisen2013; Schannwell and others, Reference Schannwell, Drews, Ehlers, Eisen, Mayer and Gillet-Chaulet2019). Moreover, because of the proximity of the Neumayer station III, numerous observations are available, e.g. ice thickness, surface velocities and most importantly surface accumulation rates, which we will use later for validation.
Figure 5. Overview of the Ekström Ice Shelf. (a) Satellite view of Ekström Ice Shelf along with location of the radar transect along the central flow line (red line) and the Kottas traverse (blue line). An independent estimate of surface accumulation via stake arrays is available on Kottas traverse, which we use to validate our results. In our model, we use the velocity data from ITS_LIVE (Gardner and others, Reference Gardner2018; Reference Gardner, Fahnestock and Scambos2022). (b) Vertical cross-section view of the radar transect, along with ice surface and base take from BedMachine Antarctica (Morlighem and others, Reference Morlighem2017), starting at the grounding line (GL). Red lines indicate four picked internal reflection horizons (IRHs). (c) Zoom in on box in B. The IRHs are numbered 1–4 in order of increasing depth. This plot is shown with the radar data used to label the IRHs in Figure I1.
4.1. Data preprocessing
First, we used Antarctic Mapping Tools (Greene and others, Reference Greene, Gwyther and Blankenship2017), BedMachine Antarctica (Morlighem and others, Reference Morlighem2017) and ITS_LIVE (Gardner and others, Reference Gardner2018; Reference Gardner, Fahnestock and Scambos2022) to obtain the surface elevation s, thickness h and velocity v for EIS. In order to define the flow tube domain for EIS, we also used the itslive_flowline tool to find two flow lines which formed the side-boundaries of the domain. The other two boundaries of the domain were the grounding line, and a straight line connecting the two flow lines. The straight line was chosen to ensure that the radar transect where data were measured is wholly contained within the flow tube domain.
We preprocessed the raw ice shelf geometry and velocity data prior to evaluating the model. This ensured numerical stability of the forward model. Using the icepack package for Python (Shapero and others, Reference Shapero, Badgeley, Hoffman and Joughin2021), we first smoothed the raw thickness data by solving a regularized minimization problem. We then solved for the best-fitting velocity by fitting a fluidity parameter in an SSA model to the observed velocity and smoothed thickness. The hyperparameters used for preprocessing are given in Appendix D.
4.2. Radar measurements of internal stratigraphy
Internal stratigraphy data along the central flow line of EIS (Fig. 5a) were acquired using a ground-based ground-penetrating radar with a center frequency of 
$50\,\mathrm{MHz}$ (pulseEKKOTM from Sensors & Software) in two consecutive field seasons (2021/22 and 2022/23) with logistic support from the Neumayer III station (Wesche and others, Reference Wesche2016; Wesche and Regnery, Reference Wesche and Regnery2022). Radar processing was done with ImpDAR (Lilien and others, Reference Lilien, Hills, Driscol, Jacobel and Christianson2020) and included trace averaging to equidistant spacing (
$10\,\mathrm{m}$), bandpass filtering (with cutoff frequencies of 20 and 
$75\,\mathrm{MHz}$), and a topographic correction using the REMA surface elevation (Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019). The latter provides observations consistent with the modeling setup. The radar detects the ice–ocean interface and continuous IRHs down to ∼200 m depth (Fig. 5b and c). Four IRHs were digitized along the entire 130 km long profile using a semi-automatic maximum tracking scheme. The vertical offset off IRHs at the profile junction in the mid-shelf region between both years is much smaller than the radar system’s wavelength in ice (
$\sim\!\! 3.4\,\mathrm{m}$). Consequently, IRHs were connected without adjustments. For the travel time-to-depth conversion, we used a depth-density profile representative for ice shelves of the Dronning Maud Land Coast (Hubbard and others Reference Hubbard(2013), eqn (1)).
4.3. Inference results
We inspect the prior over the basal melt rates as a validation of our modeling choices. The implicit prior is the same as the prior defined for the surface accumulation, with the mean shifted by the total mass balance on the flow line, 
$\dot{\mathbf{m}}$. The basal melt rate is larger (up to 
$4\,\mathrm{m}\,\mathrm{a}^{-1}$) near the grounding line and gradually stabilizes in the along-flow direction to values between 0 and 
$1\,\mathrm{m}\,\mathrm{a}^{-1}$ downstream. This is in agreement with previous estimates for basal melt profiles on this particular ice shelf (Neckel and others, Reference Neckel, Drews, Rack and Steinhage2012).
We infer the surface accumulation and basal melt rates from IRH 2 in our dataset, which has an average (ice equivalent) depth of 
$30\,\mathrm{m}$ (Fig. 6). The posterior over the surface accumulation rate has uncertainty comparable to that of the prior. However, there is a shift in the overall spatial trend of the accumulation rate; particularly, there is higher surface accumulation rate at ∼20 km from the grounding line. Accumulation rate also increases steadily downstream the flow line. As in the synthetic case, the posterior predictive distribution reproduces the observed IRH with much higher fidelity and confidence than the prior predictive distribution. The average RMSE relative to the observed IRH is 
$4.6\,\mathrm{m}$ for 1000 posterior predictive simulations, compared to 
$11.8\,\mathrm{m}$ for 1000 prior predictive simulations. The posterior predictive produces an independent estimate of the unknown age of the IRH of 
$84^{+52}_{-30}$ years (meaning a median of 84 years, and 16th and 84th percentiles of 54 and 136 years, respectively).
Figure 6. Prior and posterior (predictive) for the Ekström dataset, IRH 2, of average depth 30 m. (a and c) Prior and posterior over surface accumulation and basal melt rates respectively, starting at the grounding line (GL). Solid line is the distribution mean, the shaded region represents the 5th and 95th percentiles. (b) Cross section of the ice shelf. Prior and posterior predictive distributions for the layer closest matching the observed IRH. The vertical dashed line represents the LMI boundary for this IRH. The posterior predictive reconstructs the observed IRH with higher accuracy and lower uncertainty. The posterior predictive distribution of the age of the IRH is 
$84^{+52}_{-30}$ years (meaning a median of 84 years, and 16th and 84th percentiles of 54 and 136 years, respectively). The average RMSE relative to the observed IRH is 
$4.6\,\mathrm{m}$ for the posterior predictive distribution and 
$11.8\,\mathrm{m}$ for the prior predictive distribution.
Our method can use much deeper IRHs for the inference of accumulation and basal melt rates. For IRH 4 of the observed dataset (of average depth 
$131\,\mathrm{m}$), the proportion of the IRH that is within the LMI body is smaller. This is due to the unknown boundary condition influencing the IRH elevation at much further points along the flow line. This discarding of data has visible effects on the posteriors over the mass-balance parameters (Fig. 7). These rates are now more similar to the priors for the first 60 km of the transect and only diverge at points further down the ice shelf, where the values of accumulation and basal melt rates affect the dynamics of the IRH. Regardless, the posterior predictive distribution resembles the observed IRH at higher fidelity and precision than the prior predictive. The average RMSE relative to the observed IRH is 
$10.0\,\mathrm{m}$ for the posterior predictive distribution and 
$16.4\,\mathrm{m}$ for the prior predictive distribution. The estimated age of this IRH by our method is 
$188^{+96}_{-49}$ years. The uncertainty of the age estimates reasonably increases for deeper IRHs.
Figure 7. Prior and posterior (predictive) for the Ekström dataset, IRH 4, of average depth 113 m. Same as Figure 6 for the deeper IRH. The posterior predictive distribution of the age of the IRH is 
$188^{+96}_{-49}$ years. The average root-mean-square error relative to the observed IRH is 
$10.0\,\mathrm{m}$ for the posterior predictive distribution and 
$16.4\,\mathrm{m}$ for the prior predictive distribution.
5. Discussion
5.1. Posterior mass-balance rates are consistent between IRHs
We compare the four posteriors over the surface accumulation obtained from the Ekström IRH dataset (Fig. 8). The posteriors for the shallower IRHs 1–3 all show a similar qualitative relationship: a local maximum of the accumulation at a distance of ∼20 km from the grounding line, followed by a steady increase in the accumulation downstream. The increase in accumulation at 
$\sim\!\!\! 20\,\mathrm{km}$ is even identified in the posterior for IRH 3, despite the LMI boundary being downstream of it, at ∼30 km from the grounding line. This is reasonable, as the mass-balance parameters at a given location affect the flow field downstream of this location, and consequently, the formation of isochronal layers. For IRH 4, the LMI boundary is much further downstream at 
$\sim\!\! 60\,\mathrm{km}$. Thus, the local surface accumulation maximum at 
$\sim\!\! 20\,\mathrm{km}$ is not found; however, the overall trend of increasing surface accumulation downstream is still identified. There is a corresponding trend in the basal melt rate, as the local basal melt rate still exhibits a maximum at 
$\sim\!\! 20\,\mathrm{km}$. The reason for this is unknown, but the location corresponds both with the seaward limit of the tidal flexure zone and with the confluence region of ice originating from the eastern tributary. One or both of these factors could alter the basal melt rates inferred at this location. As we will show later (Section 5.3), this local maximum also appears in independent remote-sensing estimates.
Figure 8. Ekström Ice Shelf—dependence of posterior surface accumulation rate on depth of IRH used for inference. The posteriors are compared to the shallow layer approximation (SLA) and local layer approximation (LLA) (Waddington and others, Reference Waddington, Neumann, Koutnik, Marshall and Morse2007), and an estimate of the distribution of the accumulation rate based on measurements along the Kottas traverse. See Figure G1 for yearly Kottas measurements. As the real age of the IRHs is not known, the SBI-derived median age is used for the SLA and LLA approximations. Median ages for IRH 1–4 are 
$42^{+32}_{-12}$, 
$84^{+52}_{-30}$, 
$146^{+52}_{-38}$ and 
$188^{+96}_{-49}$ years. The LMI boundary, representing where the IRH data were masked, is shown with the brown dashed lines.
The inferred posteriors also allow us to estimate the age of the IRHs. By sampling from the posterior distribution, and evaluating the forward model with the resulting mass-balance parameter samples, we obtain a distribution of isochronal layers similar to the observed IRH, with known ages. Thus, we estimate the ages of the four IRHs as 
$42^{+32}_{-12}$, 
$84^{+52}_{-30}$, 
$146^{+52}_{-38}$ and 
$188^{+96}_{-49}$ years. This is an important finding because IRH age is otherwise only accessible through ice coring. Our results, however, depend on a realistic prior for the surface accumulation and basal melt rates, as defined in Section 2.2.3. Given a miscalibrated prior, the estimated ages would not be reliable (see Fig. H1 for an example). We hypothesize that, given an independent measurement of the IRH age, our approach could further constrain the posterior distributions over the mass-balance parameters.
The consistent spatial patterns of magnitudes of accumulation rates inferred from IRHs 1–4 are supportive of EIS being in steady state over the last hundreds of years but given that steady-state is one of our model assumptions this interpretation needs to be considered with care.
5.2. Comparison to shallow and local layer approximations
To validate our approach, we compare the inferred surface accumulation rate of our experiments with estimates from other methods. First, we computed the shallow layer approximation (SLA) and local layer approximation (LLA) as described in Waddington and others Reference Waddington, Neumann, Koutnik, Marshall and Morse(2007). Given the depth and age of IRH m, the SLA and LLA approximations for the accumulation rate 
$\dot{\mathbf{a}}$ are defined as
\begin{equation}
\begin{aligned}
\dot{\mathbf{a}}_{\text{SLA}}^{m} &= \frac{1}{\mathcal{A}_{m}}(\mathbf{s} - \mathbf{e}_{m}(\mathbf{x})),\\
\dot{\mathbf{a}}_{\text{LLA}}^{m} &= -\ln \left( 1- \frac{\mathbf{s}-\mathbf{e}_{m}(\mathbf{x})}{\mathbf{h}}\right)\frac{\mathbf{h}}{\mathcal{A}_{m}}~,
\end{aligned}
\end{equation} where 
$\mathcal{A}_{m}$ is the age of IRH m. Intuitively, the SLA takes the ice thickness above layer m and divides it by the layer age, whereas LLA accounts for strain thinning assuming a linear vertical velocity profile (which is often the case for ice-shelf flow). Since the age of the observed IRHs is not known, we use the median age of the posterior predictive distribution results. As expected, we observe that both SLA and LLA closely match the SBI posterior mean accumulation rate for the shallow IRHs of median estimated ages 42 and 84 years (Fig. 8). As the strain rates of the flow are small, the relatively shallow IRHs (mean ice equivalent depth of 
$30\,\mathrm{m}$) have not notably deformed, and hence the assumptions of SLA and LLA are appropriate. However, for the deeper IRHs 3 and 4 of estimated ages 146 and 188 years, we see that both SLA and LLA estimates diverge from our posterior mean accumulation rate. This shows that more involved approaches are required when using deeper IRHs for inference. For deeper IRHs where the SLA and LLA no longer applied, Steen-Larsen and others Reference Steen-Larsen, Waddington and Koutnik(2010) inferred the surface accumulation rates on grounded ice using a Monte Carlo approach. By treating the age of the IRH as an additional parameter to infer, they were able to identify the age of the IRH with high confidence. Extensions of our approach could incorporate this parameterization to reduce the uncertainty of the inferred IRH.
5.3. Comparison with independent estimates of surface accumulation and basal melting
For the Ekström transect comparison data are provided by repeat readings of accumulation stakes in 
$500\,\mathrm{m}$ spacing along the nearby Kottas traverse (Fig. 8). Yearly readings are available in the period 1996–2005 and on a yearly to three-yearly interval between 2014 and 2023 (Mengert, Reference Mengert2018). We use this dataset to construct a direct estimate of time-averaged surface accumulation rate along the central flow line transect over these periods. For this, we project the measurements from the Kottas traverse to the flow line transect, taking into account an increased uncertainty for increasing projection distance (see Appendix G for details).
The Kottas traverse accumulation measurements closely match the posterior means of our approach (Fig. 8) for IRHs 1 and 2. As the accumulation rate measurements on the Kottas traverse span the past 26 years, it may not be a good validation for the deeper IRHs. Regardless, the Kottas accumulation rate measurements lie within the posterior uncertainty for IRHs 3 and 4. These comparisons further corroborate our approach and highlight the advantages of uncertainty-aware methods, especially as the measured accumulation rates also varied considerably year-to-year (Appendix G).
We also compare our inferred basal melt rates with independent measurements of basal melt rates. In Fig. 9, we show the basal melt rates inferred from IRH 1 for the Ekström dataset, in comparison to independent estimates of basal melt rates through satellite altimetry data (Adusumilli and others, Reference Adusumilli, Fricker, Medley, Padman and Siegfried2020), and through airborne radar measurements of the ice-shelf thickness (Neckel and others, Reference Neckel, Drews, Rack and Steinhage2012). We observe a quantitative match between our posterior basal melting rate and the estimates from Adusumilli and others Reference Adusumilli, Fricker, Medley, Padman and Siegfried(2020). The ice-shelf wide melt rate estimates show overall comparatively little spatial variability. The most notable difference is that our basal melt rates show more variability in the first 20 km of the profile, and this could be due to the proximity of the grounding zone where the SSA approximation does not hold. However, we note that also the satellite-derived estimates show this oscillation in basal melt rates albeit with a smaller magnitude. The estimates from Neckel and others Reference Neckel, Drews, Rack and Steinhage(2012) excluded the grounding zone area but otherwise show a good match in magnitude but with much less spatial variability. This is because they decided to apply spatial smoothing the degree of which could be revisited given the new results derived here.
Figure 9. Basal melting rates comparison. (a) Map of basal melt rates for Ekström Ice Shelf, using data from Adusumilli and others Reference Adusumilli, Fricker, Medley, Padman and Siegfried(2020). (b) Comparison of inferred basal melt rates from IRH 1 to independent estimates of basal melt rates, calculated on the flow line transect.
The good quantitative match with independently collected data both for surface accumulation and basal melting increases our confidence for our inferred surface accumulation and basal melt rates. However, these results are limited by some of our modeling assumptions, which we discuss next.
5.4. Limitations of modeling approach
The fidelity at which the posterior predictive distributions reproduce the observed IRHs of EIS (Figs. 6 and 7) supports our modeling choices for this ice shelf, as the combination of the forward model and accumulation rate prior distribution are sufficiently expressive to reproduce the IRHs.
However, our inferred surface accumulation and basal melt rates rely on the modeling assumption that the internal radar data are collected on a flow-line transect. This assumption is required to conclude that the ice observed in the internal stratigraphy is indeed the accumulated ice modeled in our domain, as commonly assumed in the literature (Waddington and others, Reference Waddington, Neumann, Koutnik, Marshall and Morse2007; Steen-Larsen and others, Reference Steen-Larsen, Waddington and Koutnik2010; Theofilopoulos and Born, Reference Theofilopoulos and Born2023). Similarly, in the context of ice rises, it is often assumed that the transverse velocity is negligible relative to the vertical velocity of the ice, so that the surface accumulation rate can be inferred along the same transect as the IRH data (Callens and others, Reference Callens, Drews, Witrant, Philipp and Pattyn2016; Koch and others, Reference Koch2024). However, many of the available IRH measurements do not align with flow lines of the ice sheet. In this case, our assumption would not be valid.
Because our observations are on a flow line transect only, it is difficult to judge to what extent unidentified three-dimensional (3-D) effects overprint our analysis. Previous approaches (Pattyn and others, Reference Pattyn2012) have had similar limitations because ice-flow divergence and/or convergence could not be predicted by their 2-D forward model. They concluded that their inferred basal melt rates which best matched the radar stratigraphy would be a lower boundary because ice flow on the ice shelf was convergent. In our case, we do correct for the observed convergence from observed surface velocities along the flow tube. The normal component of ice-flow is always 
$ \lt 20\,\mathrm{m}\,\mathrm{a}^{-1}$ and often 
$ \lt 1\,\mathrm{m}\,\mathrm{a}^{-1}$. This is small compared to the along-flow velocities and also compared to the total surface mass balance accumulated along the flow tube. Together with the empirical validation with independently collected surface accumulation and basal melt rates, this increases our confidence that our modeling approach yields trustworthy results. Yet, a more rigorous quantification of 3-D effects, for example, in a synthetic study using a 3-D forward model can provide further validation. Others have made progress in this direction, for example, Wolovick and others Reference Wolovick, Moore and Zhao(2021) consider a 3-D steady-state ice sheet and jointly infer the temporally averaged accumulation rate and geothermal heat flow. This is done by using more radar attributes in addition to stratigraphy such as existence or absence of subglacial water and/or basal freeze-on.
Common to most previous studies is the steady-state assumption which is imposed because the inclusion of transient ice thickness changes increases the model parameter space to a degree which cannot easily be solved in the inverse problem, particularly with quantified uncertainties over the inferred parameters. Ways forward in this regard could be deterministic gradient descent schemes with explicitly calculating sensitivity matrices as suggested by Theofilopoulos (Reference Theofilopoulos2022); Theofilopoulos and Born (Reference Theofilopoulos and Born2023) and this will be an important step forward to better exploit the growing IRH archive for ice-sheet modeling.
5.5. SBI as a tool for geoscientific inversion problems
The inverse problem tackled in this work typifies geoscientific inverse problems, as the forward model is defined in terms of a partial differential equation, and the parameters are high dimensional and vary in space. Hence, it is valuable to compare the SBI approach in this case to the wide variety of methods and algorithms that have been developed to solve geoscientific inverse problems. In the remainder of the section, we discuss NPE as used in our work. However, there exist other variants of SBI with relative advantages and disadvantages, depending on the problem setting (Cranmer and others, Reference Cranmer, Brehmer and Louppe2020; Lueckmann and others, Reference Lueckmann, Boelts, Greenberg, Goncalves and Macke2021). In particular, we provide a brief discussion of the neural likelihood estimation (NLE) variant in Appendix C.4.
The SBI approach as presented here has two key features. First, we estimate the Bayesian posterior distribution, providing quantitative uncertainty estimates. Modeling uncertainty is important as it can highly influence and propagate to future modeling predictions. Additionally, locations of high uncertainty show areas requiring further study, helping to guide future work. This is in contrast to deterministic inversion methods, which do not estimate uncertainty, or likelihood-based inference methods which are not possible when the likelihood defined by the simulator is not known. Thus, approximate Bayesian methods and SBI in particular can be applied to a larger class of inference problems. Second, a unique advantage of single round SBI methods (Cranmer and others, Reference Cranmer, Brehmer and Louppe2020) such as NPE (as used in our study) is amortization. Our method as presented here is not yet fully amortized, instead amortizing the vast majority of the computational cost, as preprocessing relies on the observed value of X. In order to train the density estimator 
$q_{\phi}(\dot{\mathbf{a}}|X_{k}^{m})$, we first calculate 
$X_{k}^{m}$ for each simulation dependent on the value of 
$X_{o}^{m}$. Our method still amortizes the cost of simulating the forward model many times, which is by far the largest computational cost in the approach. In the Ekström example, we have evaluated the forward model a total of 190 000 times, accounting for ∼99% of the total computation cost (Appendix E). This amortization is specific to the geometry and velocity of the EIS; different geometries and velocities change the dependence of the internal layers on the mass-balance parameters, which would require simulating from a new model.
On the other hand, SBI faces some limitations as an inference tool. Primarily, SBI methods are known to require a large number of simulations to be trained (Lueckmann and others, Reference Lueckmann, Boelts, Greenberg, Goncalves and Macke2021). This problem suffers from the curse of dimensionality—the number of simulations required scales exponentially with the number of parameters we are trying to infer (in this work, we limited the number of parameters to 50). This is particularly challenging for geoscientific problems, where typically the parameters of interest vary spatially (and temporally), and thus the number of parameters can grow very large. The SBI approach needs to be adapted to more efficiently represent high-dimensional, spatially varying parameters θ at high resolutions. Some potential approaches are polynomial or spectral representations. Future work should also explore variants of SBI that are better suited to high-dimensional or even continuous parameters (Ramesh and others, Reference Ramesh2022; Geffner and others, Reference Geffner, Papamakarios and Mnih2023). Finally, SBI works under the assumption that the forward model is well-specified, meaning that given samples from the prior, it can generate simulations closely resembling the observation. The posteriors obtained by SBI can be strongly biased when this is not the case (Cannon and others, Reference Cannon, Ward and Schmon2022). Work to address this concern has been done, e.g. by incorporating the model mismatch into the forward model (Ward and others, Reference Ward, Cannon, Beaumont, Fasiolo and Schmon2022), as done in our work using the calibrated noise model. However, designing and calibrating such noise models for each inference task are challenging, and a standard approach for addressing model mismatch does not yet exist.
6. Conclusions
We presented a novel approach for inferring the spatially varying surface accumulation and basal melt rates along ice-shelf flow lines from radar measurements of their internal stratigraphy. We validated the method on a synthetic ice shelf example and inferred the surface accumulation and basal melt rates along a flow line in EIS, Antarctica. We separately inferred the mass-balance parameters from four different IRHs. The inferred distributions were further validated by independent stake array measurements of surface accumulation rates uniquely available in Ekström Ice Shelf. Using our approach, we were able to estimate the otherwise unknown age of the IRHs as 
$42^{+32}_{-12}$, 
$84^{+52}_{-30}$, 
$146^{+52}_{-38}$ and 
$188^{+96}_{-49}$ years. The presented approach can be transferred to other Antarctic ice shelves and also to other flow regimes such as grounded ice. A strength of our approach is the principled uncertainty estimates in the inferred surface accumulation and basal melt rates. These uncertainty estimates can be integrated in future projections of the Antarctic ice sheet (Verjans and others, Reference Verjans, Robel, Seroussi, Ultee and Thompson2022; Ultee and others, Reference Ultee, Robel and Castruccio2024). We identified avenues for future work as more can be learned by relaxing the steady-state assumption on the ice shelf. The forward model and inference framework should be adapted to account for potential transient signals in the mass-balance parameters.
This work was an example use case of SBI for a geoscientific inverse problem. We showcased the strengths of SBI as a likelihood-free approach to approximate the Bayesian posterior, amortizing the cost of simulating the forward model many times. SBI can become more applicable to such inverse problems involving spatially (and temporally) varying parameters if it can be extended to deal with the challenge of high-dimensional parameter inference.
Finally, our approach highlights the value of internal stratigraphy measurements. Initiatives to map the Antarctic-wide internal stratigraphy (e.g. Bingham and others, Reference Bingham2024) can provide invaluable data toward uncovering the history of the Antarctic ice sheet. Sophisticated inference methods could be combined with such a dataset to provide a new, independent, Antarctica-wide parameterization of accumulation and basal melt rate histories.
Data availability statement
The extracted IRH elevations along the Ekström transect are available in Oraschewski and others Reference Oraschewski, Moss, Koch, Ershadi, Eisen and Drews(2024a). The processed radar data are additionally available in Oraschewski and others Reference Oraschewski, Moss, Koch, Ershadi, Eisen and Drews(2024b). Simulation data available in Moss and others Reference Moss(2023).
Software availability
Code for preprocessing Ekström Ice Shelf data and generating synthetic ice shelf data is available in Moss and others Reference Moss(2024a). Code for layer tracing forward model and simulation-based inference workflow is available in Moss and others Reference Moss(2024b).
Acknowledgements
The authors would like to thank Daniel Shapero for his inputs on use of icepack. The authors would also like to thank Andreas Born and Therese Riekch for insightful discussions on the implementation of the layer tracing solver for calculating internal stratigraphy. We acknowledge excellent logistic support from staff at Neumayer Station III and the GrouZe team on-site.
This work was funded by the German Research Foundation (DFG) under Germany’s Excellence Strategy—EXC number 2064/1—390727645 and SFB 1233 ‘Robust Vision’ (276693517) and the German Federal Ministry of Education and Research (BMBF): Tübingen AI Center, FKZ: 01IS18039A. Reinhard Drews and Vjeran Višnjević were supported by an Emmy Noether grant of the Deutsche Forschungsgemeinschaft (DR 822/3-1). We acknowledge the support by the German Academic Scholarship Foundation to Falk M. Oraschewski. Field observations were supported by the Alfred Wegener Institute through logistic grants AWI_ANT_23 (Drews) and AWI_ANT_8 (Eisen). We acknowledge support from the Open Access Publication Fund of the University of Tübingen. Guy Moss is a member of the International Max Planck Research School for Intelligent Systems (IMPRS-IS).
Appendix A. Forward model details
The layer tracing scheme described in Section 2.1 is equivalent to solving a set of advection equations for a set of layers. Here, we explicitly write down the advection equations solved and the boundary conditions defined. We correct for 2-D effects in the advection equations by adding a term for the normal flow into the flow line. We then account for the inflow boundary condition at the grounding line x = 0.
A 1-D advection equation for a layer on a flow line reads
\begin{equation}
\frac{\partial h_{l}}{\partial t} = v_{x}\frac{\partial h_{l}}{\partial x}.
\end{equation}In practice, real flow lines have some incoming or outgoing (normal) flux, qy, where y denotes the horizontal direction perpendicular to x. We account for this normal ice flux and instead solve
\begin{equation}
\frac{\partial h_{l}}{\partial t} = v_{x}\frac{\partial h_{l}}{\partial x} + r_{l}(x)\frac{\partial q_{y}(x)}{\partial y},
\end{equation}where 
$r_{l}(x) = h_{l}(x)/h(x)$ is the ratio of thickness of layer l to the total thickness of the shelf. This equation holds due to the plug flow assumption, in which the flux divergence 
$\partial q_{y}/\partial y$ is independent of the depth z. The quantity 
$\partial q_{y}/\partial y$ is constant for all layers and independent of the layer thickness. This normal flux component accounts for lateral compression or extension of the flowtube centered on the flow line, and in our case we estimate it from satellite inferred velocities. For EIS, the normal flow component is small compared to the total mass balance along the flow lines, but for other cases this correction might be much more significant.
In order to define the inflow boundary condition, we would need to know the relative thickness of the incoming layers (or alternatively, the vertical age distribution at x 0). This is typically not available in radar measurements of the stratigraphy. It is, therefore, important to use only the IRH elevation information within the LMI body of the flow line. This is the region of our domain which is independent of the boundary condition we chose at x 0 (Fig. A1). This region can be found by tracking the trajectory traced by a particle initially at 
$(x=x_{0},z=s(x_{0}))$. The LMI body is the region of the domain above this path. As a consequence of this consideration, we need to discard more of the IRH elevation data for the deeper IRHs in the dataset.
Figure A1. Local meteoric ice (LMI) body. The layer elevations in the LMI body (shaded region) are independent of the inflow boundary conditions. Outside the LMI body, the layer elevations are dependent on this boundary condition, and so using IRH observations within this region would require assuming the internal stratigraphy of the incoming ice.
For a complete definition of the simulator, we still need to define some boundary condition at x 0. Since the true boundary condition is not known, we choose an inflow boundary condition which improves the numerical stability of the layer tracing model,
\begin{equation}
\frac{\partial h_{i}}{\partial t} = \left.\frac{\partial q_{x}}{\partial x}\right\vert_{x=0} + r_{i}(x)\frac{\partial q_{y}(x)}{\partial y}.
\end{equation}This boundary condition makes the implicit assumption that near the inflow boundary, the layer thickness profile is a scalar multiple of the total thickness,
\begin{equation}
h(x)\frac{\partial h_{i}(x)}{\partial x} = h_{i}(x)\frac{\partial h(x)}{\partial x}.
\end{equation}We define a similar boundary condition at 
$x=L_{x}$; however, this boundary condition does not have a similar effect on the applicability of the IRH data.
Appendix B. Calibrating the simulator
We calibrate hyperparameters of the noise model using a set of 1000 simulations for the same ice shelf geometry, velocity and mass-balance parameters prior. This set of calibration simulations is not used again to train NPE (Section 2.2.1) to avoid overfitting. We use the same set of calibration simulations to calibrate the per-layer LMI boundary mask (Section 2.2.2). Algorithm 1 defines our calibration procedure for the parameters of the noise distribution for EIS.
Algorithm 1: Noise model calibration
In Section 2.2.3, we used the LMI boundary i(m) for each IRH m in order to define the domain 
$\mathbf{x}_{i\geqslant i(m)}$ on which we compare the observed IRH data to the simulated isochronal layers. The physical interpretation of the LMI boundary is found in Appendix A, and here we specify the exact definition of i(m).
For each simulation k in the calibration set, we define the trajectory of a particle starting at the surface of the inflow boundary 
$(x=0,z=s(0))$ by 
$p^{k}(t) = (x^{k}(t),z^{k}(t))$. In practice, since we are on a flow line, 
$v_x \gt 0$ everywhere, and thus x k increases monotonically with t. Therefore, the trajectory traces out a unique curve 
$z^{k}(x)$ in the domain.
Using this path, we define the LMI boundary 
$i(m,k)$ for simulation k and IRH m to be the first point in the x domain such that the path is below the IRH elevation,
\begin{equation}
i(m,k) = \min_{i} \{i : z^{k}(x_{i}) \lt e_{m}(x_{i}) \}.
\end{equation}In the case that the path stays above the IRH for the entire domain, we define 
$i(m,k) = L_{x}$, meaning that the IRH is entirely outside the LMI body.
This definition gives a different LMI body for each simulation. In order to perform inference, we require one a fixed LMI body boundary to use across all simulations. We define i(m) ‘pessimistically’ as the 
$75{\text{th}}$ percentile of the calculated 
$i(m,k)$ boundaries in the calibration set. Therefore, some of the simulated layer elevations we use will be dependent on the unknown boundary conditions but a small amount that should not affect the inferred results.
Appendix C. Additional information on simulation-based inference
C.1. Tractable and intractable likelihood functions
Here, we provide intuition about the difference of forward models which have intractable and tractable likelihood functions. Forward models which have tractable likelihood functions are ones where the density 
$p(X|\theta)$ can be explicitly evaluated for a given 
$(X,\theta)$ pair, where X refers to the observed data and θ are the parameters of the model. As an example, a common setting in glaciology is a deterministic forward model with additive observational noise. Such forward models are of the form
\begin{equation}
X = f(\theta) + \epsilon ,
\end{equation} where 
$f(\theta)$ is a deterministic function, and ϵ follows some distribution, 
$\epsilon\sim p_{\epsilon(\epsilon)}$. In such models, the likelihood 
$p(X|\theta)$ can be evaluated as
\begin{equation}
p(X|\theta) = p_\epsilon (X-f(\theta)).
\end{equation}A more general setting considers nondeterministic forward models, where latent variables z are sampled, and then affect further computations in the forward model. As a motivational example, consider a stochastic differential equation solved with an Euler–Maruyama scheme, that is, an iterative update
\begin{equation}
X_{t+1} = X_t + f(X_t,\theta, t)dt + g(X_t,\theta,t)z_t\sqrt{dt},
\end{equation} where f is now a deterministic function dependent on the state Xt and parameters θ and t, g is a deterministic function, and 
$z_t \sim \mathcal{N}(0,I)$ is random noise. Suppose the outcome from this model is the state at the end of the simulation, 
$X=X_T$. In order to compute the likelihood of this model, we need to integrate over all the noise samples 
$z = [z_1,\ldots,z_{T-1}]$,
\begin{equation}
p(X|\theta) = \int p(X,z|\theta)dz.
\end{equation} In general, approximating this integral is significantly more expensive than generating a sample from 
$p(X|\theta)$ by simulating, and so we call this likelihood intractable
For our forward model, the intractability of the likelihood follows from the fact that we select the closest layer to the IRH after sampling the noise. Therefore, the noise model cannot be efficiently calculated as in Eqn (C2). Note that, even in the case of tractable likelihood functions, it might be beneficial to use ‘likelihood-free’ inference methods, for example in the case that the forward model is computationally expensive to compute.
C.2. Normalizing flows
While NPE can be performed with any valid parameterization of the variational distribution 
$q_\phi(\theta|X)$, in this work, we make the common choice to use a normalizing flow model. In the following, we adapt the notation of Papamakarios and others Reference Papamakarios, Nalisnick, Rezende, Mohamed and Lakshminarayanan(2019a).
The goal of normalizing flows is to learn a map from a distribution 
$u\sim p_u(u)$ to a distribution of the same dimensionality 
$x\sim p_x(x)$. Here, 
$p_x(x)$ is the target distribution, and 
$p_u(u)$ is a simple distribution that can easily be evaluated and sampled from (e.g. a multivariate Gaussian). If we can express x as a differentiable, invertible transformation of u, i.e. 
$x = T(u)$, then we can also evaluate the probability
\begin{equation}
p_x(x) = p_u(u)|\text{det} J_T(u)|^{-1},
\end{equation} where 
$u = T^{-1}(x)$ and JT is the Jacobian matrix of the transformation T.
Normalizing flows are then defined by a sequence of transformations 
$T_1,\ldots,T_K$ which transform 
$z_0=u$ to 
$z_K = T_K\circ \cdots \circ T_1 (u) = x$. Each of the transformations is learnable. The composition of differentiable, invertible transformations is also invertible and differentiable, with the determinant of the Jacobian satisfying
\begin{equation}
\text{det} J_{T_2\circ T_1}(u) = \text{det} J_{T_2}(T_1(u)) \cdot \text{det}J_{T_1}(u).
\end{equation} One common way to parameterize the individual transformations Ti is through coupling transforms, where the input z is split into two parts, 
$z = [z_{1:d},z_{d+1:D}]$, where D is the dimensionality of 
$z\in\mathbb{R}^{D}$. The first part, 
$z_{1:d}$, remains unchanged. The second part is transformed elementwise, 
$z'_{i} = \tau(z_i;h_i)$, where τ is some monotonic function of zi conditioned on 
$h_i = F(z_{1:d})$, for some learnable function F of the first d components of z. The vector is permuted between each layer, so that not the same components are mapped through the identity with each transformation. An advantage of this parameterization is that the Jacobian matrix can be calculated as the product of its diagonal components, making the evaluation of the log-probabilities significantly faster. In this work, we use neural spline flows (NSFs, Durkan and others Reference Durkan, Bekasov, Murray and Papamakarios2019). In NSFs, τ are monotonic spline functions, which are analytically invertible, yet highly flexible.
C.3. Connection between NPE loss and Kullback-Leibler divergence
We note that the expected forward KL divergence between the true posterior 
$p(\theta|X)$ and the variational distribution 
$q_\phi(\theta|X)$ can be decomposed as
\begin{equation}
\begin{aligned}
&{\mathbb{E}}_{p(x)}[D_{\text{KL}}(p(\theta|X||q_\phi(\theta|X)))] = \mathbb{E}_{p(X)p(\theta|X)}\left[\log \frac{p(\theta|X)}{q_\phi(\theta|X)}\right]\\
&= -\mathbb{E}_{p(X)p(\theta|X)}[\log q_\phi(\theta|X)] + \mathbb{E}_{p(X)p(\theta|X)}[\log p(\theta|X)],
\end{aligned}
\end{equation}where the first term on the right hand side is the negative of the loss in Eqn (8), and the second term is independent of the variational parameters ϕ and is thus a constant with respect to the variational parameters ϕ which we optimize. Therefore, minimizing the NPE loss (Eqn (8)) is equivalent to minimizing the expected forward KL divergence between the variational and true posterior distributions.
C.4. Neural Likelihood Estimation
In our work, we use NPE to directly approximate the posterior distribution. However, other variants of SBI exist. In particular, a prominent method of SBI is NLE (Wood, Reference Wood2010; Papamakarios and others, Reference Papamakarios, Sterratt and Murray2019b). In NLE, the goal is to approximate the likelihood function 
$p(X|\theta)$ from a dataset of simulations 
$\{\theta_k, X_k\}_{k=1}^{K}$. This likelihood can then be used to perform Bayesian inference using the existing set of Bayesian inference methods, such as Markov chain Monte Carlo (MCMC), or variational inference (Blei and others, Reference Blei, Kucukelbir and McAuliffe2017).
Similarly to NPE, we define a variational distribution (e.g. a normalizing flow), 
$q_\phi(X|\theta)$, which can be trained by minimizing the negative log-probabilities of the model outcomes Xk given model parameters θk
\begin{equation}
\mathcal{L}(\phi) = \mathbb{E}_{\theta_k\sim p(\theta),X_k\sim p(X|\theta_k)}[ -\log q_{\phi}(X_k|\theta_k)].
\end{equation} Notice that the conditioning of 
$\theta_k,X_k$ is reversed relative to Eqn (8).
The likelihood model learned in NLE can also be thought of as an emulator of the forward model, as it allows us to draw samples 
$X\sim p(X|\theta)$. This makes NLE methods advantageous also for problems where the likelihood is tractable, as the learned emulators can typically be evaluated faster than the original forward model. Related ideas have been explored for glaciological problems. For example, Tarasov and others Reference Tarasov, Dyke, Neal and Peltier(2012) use Bayesian neural networks as a surrogate model for particular parameter-output pairs in a glacial systems model, which they then use to perform inference with MCMC. Similarly, Brinkerhoff and others Reference Brinkerhoff, Aschwanden and Fahnestock(2021) use a deterministic residual neural network as a surrogate model to predict a low-dimensional representation of outcomes from a subglacial hydrology model. They then derive an approximation of the uncertainty in this model to arrive at a nondeterministic likelihood which can then be used for inference.
Appendix D. Further implementation details
For completeness, we give the values of the important hyperparameters involved in the workflow. We provide the details of the spatial and temporal resolutions of our various simulations, along with the regularization strengths of the smoothing of the EIS data (Tables D1, D2 and D3).
Table D1. Hyperparameters for synthetic ice shelf spin-up modeling
Table D2. Hyperparameters for the preprocessing of the data for Ekström Ice Shelf
Table D3. Layer tracer forward model simulation configuration
During inference, we applied NPE as implemented in the sbi package (Tejero-Cantero and others, Reference Tejero-Cantero2020) to obtain the results for both the synthetic (Section 3) and Ekström (Section 4) ice shelves. We used the NSF (Durkan and others, Reference Durkan, Bekasov, Murray and Papamakarios2019) as implemented in Lueckmann and others Reference Lueckmann, Boelts, Greenberg, Goncalves and Macke(2021), and with the same architecture of five flow transformations, two residual blocks of 50 hidden units each, ReLU nonlinearity and 10 bins. We also embedded the 500-dimensional observation of layer elevations to a 50-dimensional summary statistic used as the condition for the NSF. The embedding network consisted of two convolutional layers with kernel size 5, each followed by ReLU activations and max pooling with kernel size 2. The number of output channels for the two convolutional layers were 6 and 12, respectively. The output channels of the second convolutional layer were then concatenated and fed through two fully connected linear layers, each followed by ReLU activations. The number of hidden units was set to 50. Training was done as in Lueckmann and others Reference Lueckmann, Boelts, Greenberg, Goncalves and Macke(2021), with the exception that the batch size was set to 1000 (default is 50). For each NPE run, we train five networks initialized with different random seeds and report in our results the run with the best validation loss.
Following the work of Lueckmann and others Reference Lueckmann, Boelts, Greenberg, Goncalves and Macke(2021), we split the simulation dataset into training and validation datasets with a 90–10 split. We optimize the loss in Eqn (8) using an Adam Optimizer (Kingma and Ba, Reference Kingma and Ba2015) with a batch size of 50, a learning rate of 0.0005, with the maximum gradient norm clipped to 5.0. For each training epoch, we calculate the validation loss on the entire validation dataset. If the validation loss has not surpassed its best value for 20 training epochs, we assume convergence and stop training.
Appendix E. Computational costs
We provide a breakdown of the approximate computational costs of the different stages in our workflow for both synthetic and Ekström Ice Shelves in Tables E1 and E2, respectively. These are dependent the hardware used and vary stochastically as a result of random number generators. This section is intended to provide intuition into the relative scales of the different stages of the workflow, rather than exact measurements. We had access to 16-core Intel Xeon Gold 2.9 GHz CPU nodes and Nvidia RTX 2080ti GPU nodes. While the large number evaluations of the forward layer tracing model were by far the most computationally intensive section of the workflow in both cases, these simulations were trivially performed in parallel across CPU cores, thus reducing the wall-clock time of the workflow.
Table E1. Synthetic ice shelf approximate computational cost breakdown. Some tasks are embarrassingly parallelizable—parallel resources and times are shown in square brackets. All times reported in minutes
Note: Only the last two tasks need to be repeated for each IRH measurement.
Table E2. Ekström Ice Shelf approximate computational cost breakdown. Some tasks are embarrassingly parallelizable—parallel resources and times are shown in square brackets. All times reported in minutes
Note: Only the last two tasks need to be repeated for each IRH measurement.
This analysis highlights the advantages of our amortized approached to inference. For EIS, the total computational time of the noiseless simulations accounted for 
$\approx 99.8\%$ of the total computational time of the workflow. The simulations need not be repeated when we infer from other IRHs, greatly benefiting the computational efficiency of inference as the number of IRHs in the dataset increases. This advantage is slightly reduced when considering the parallelization we have used (Table E2), as the training of the probability density estimator is not easily parallelizable across GPU nodes. Accounting for this still results in 
$\approx 97.6\%$ of the computational cost being amortized.
Appendix F. Additional results
We show the posterior and posterior predictive distributions when inferring from isochronal layer 3 in the synthetic ice shelf dataset, of age 150 years (Fig. F1). This isochronal layer has average depth of 
$120\,\mathrm{m}$, comparable to the deepest IRH in the Ekström dataset. While the uncertainty is much higher than for the shallow layer, the posterior over the surface accumulation still shows a higher mean than the prior in the downstream section of the ice shelf. Additionally, the posterior predictive reconstructs the ground truth isochronal layer better than the prior predictive. The mean RMSE of the posterior predictive layers relative to the ground truth is 
$13.6\,\mathrm{m}$, compared to 
$19.8\,\mathrm{m}$ for the prior. Therefore, we are still able to reconstruct additional information about the mass balance parameters, even from much deeper layers.
Figure F1. Prior and posterior (predictive) for the synthetic dataset. (a and c) Prior and posterior over surface accumulation and basal melt rates respectively for layer 3 of the synthetic ice shelf, of age 150 years. Solid line is the distribution mean, the shaded region represents the 5th and 95th percentiles. The ground truth (GT) parameters used to generate the reference isochronoal layer are also shown. (b) Cross section of the ice shelf. Prior and posterior predictive distributions for the layer closest matching the ground truth isochronal layer. The vertical dashed line represents the LMI boundary for this isochronal layer. The posterior predictive reconstructs the observed layer with higher accuracy and lower uncertainty. The posterior predictive distribution of the age of the isochronal layer is 
$224^{+139}_{-71}$ years.
We also explore the dependence of the inferred posterior over surface accumulation rate on the depth of the layer used for inference in the synthetic case (Fig. F2), similar to the analysis done in Section 5.2 for EIS. For the synthetic ice shelf, the surface accumulation and basal melt rates were held constant for the entire simulation time, and hence the increased uncertainty with depth seen in Fig. F2 highlights that information about the mass balance parameters is gradually lost with time as a result of the action of the simulator. Indeed, for the deepest ground truth isochronal layer of average depth 
$183\,\mathrm{m}$, the posterior distribution is almost identical to the prior distribution.
Figure F2. Synthetic shelf: dependence of posterior surface accumulation rate on depth of layer used for inference.
Finally, we report the RMSE in the predicted isochronal layer elevations, relative to the true IRH (for EIS) or the ground truth isochronal layer (for the synthetic ice shelf). This is done across simulations from 1000 simulations for each of the prior and posterior distributions, for each IRH. The RMSE is consistently lower for the posterior predictive distribution for all depths, for both the synthetic ice shelf (Table F1) and EIS (Table F2).
Table F1. Synthetic ice shelf—prior and posterior predictive distribution root-mean-square error (RMSE) relative to ground truth IRH, estimated from 1000 samples. The mean and standard deviations (SDs) in the RMSE are reported. All values are in meters
Table F2. Ekström Ice Shelf—prior and posterior predictive distribution root-mean-square error (RMSE) relative to ground truth IRH, estimated from 1000 samples. The mean and standard deviations (SDs) in the RMSE are reported. All values are in meters
Appendix G. Kottas traverse data
Here we describe the mapping of the surface accumulation measurements on Kottas traverse to the flow line transect. For each measurement year, and each location 
$\tilde{x}_{i}$ on the Kottas traverse, we find the nearest point xi on the flow line transect. We assume the accumulation rate at this point to be normally distributed, with mean 
$\tilde{\dot{a}}_{i}$ (the Kottas traverse measurement at 
$\tilde{x}_{i}$), and variance 
$\sigma_{\dot{a}}^{2}||\tilde{x}_{i}-x_{i}||^{2}_{2}/l_{\dot{a}} $. We set the length scale 
$l_{\dot{a}} = 2.5\,\mathrm{km}$ and the accumulation rate variance 
$\sigma_{\dot{a}}^{2} = 0.25^{2}\,\mathrm{m}^{2}\,\mathrm{m}^{-2}$. These values are chosen in accordance to the definition of the surface accumulation rate prior distribution (Section 2.2.3).
The yearly variations of the estimated surface accumulation as measured using the stake line along the Kottas traverse (Fig. G1) reasonably agree with the posterior inferred using IRH 2. These show that there is a high year-to-year variability in the surface accumulation (Mengert Reference Mengert2018, Fig. 13) even in this steady-state region in East Antarctica. However, the mean of these yearly measurements matches the inferred posterior mean closely.
Figure G1. Yearly variations of the Kottas surface accumulation stakes measurement dataset. Years shown are 1995–2005 and 2017–19. These are compared with the posterior distribution inferred using IRH 1 of the Ekström IRH dataset.
Appendix H. Synthetic results with miscalibrated prior
An additional outcome of our approach is the estimation of the age of isochronal layers. However, the validity of this estimate depends on the prior distribution containing the true mass-balance parameters. When the true mass-balance parameters have low probability under the prior distribution, the resulting estimate for the age of the isochronal layer can be wrong. We test this on the synthetic ice shelf by choosing a ground truth surface accumulation 
$\dot{\mathbf{a}}_{\text{mis}}$ that has low probability under the prior distribution. We calculate the isochronal layer of age 100 years under this ground truth and obtain the posterior distribution as before, using the same set of simulations as in the main text. The posterior does not capture that the mean surface accumulation rate should be higher than what is defined in the prior (Fig. H1a). However, this is not a failure of the inference method, as we can see that the posterior predictive still reconstructs the ground truth isochronal layer at higher fidelity than the prior (Fig. H1b). The predicted age of the isochronal layer 
$164^{+101}_{-44}$ years, which greatly overestimates the true age of 100 years.
Figure H1. Prior and posterior (predictive) for the synthetic ice shelf with the low-probability ground truth. (a and c) Prior and posterior over surface accumulation and basal melt rates respectively for an isochronal layer of age 100 years. Solid line is the distribution mean, the shaded region represents the 5th and 95th percentiles. The ground truth (GT) parameters used to generate the reference isochronal layer are also shown. (b) Cross section of the ice shelf. Prior and posterior predictive distributions for the layer closest matching the ground truth isochronal layer. The vertical dashed line represents the LMI boundary for this isochronal layer. The posterior predictive reconstructs the observed layer with higher accuracy and lower uncertainty. The posterior predictive distribution of the age of the isochronal layer is 
$164^{+101}_{-44}$ years.
Appendix I. Radar data
Figure I1. Radargram along transect. Zoom in on section of vertical cross-section view of the radar transect (Figure 5c). Color gradient indicates radargram data from radar survey of transect. The four labeled IRHs picked from this radargram are labeled in order of depth.
Open access
