Ice dynamics

Viscous flow

The flow of the ice sheet over a volume $\Omega$ is modeled as a low Reynolds number fluid using a hydrostatic approximation to Stokes’ equations (Pattyn, Reference Pattyn2003)

(1)$$\nabla \cdot \tau' = \rho_i g \nabla z_{\rm s}\comma\; $$

where

(2)$$\tau' = \left[\matrix{2\tau_{xx} + \tau_{yy} & \tau_{xy} & \tau_{xz} \cr \tau_{xy} & \tau_{xx} + 2\tau_{yy} & \tau_{yz} }\right].$$

$z_{\rm s}$ is the glacier surface elevation, $\rho _i$ is the ice density, $g$ is the gravitational acceleration and $\tau _{ij}$ is a component of the deviatoric stress tensor given by

(3)$$\tau_{ij} = 2\eta \dot{\epsilon}_{ij}\comma\; $$

with $\dot {\epsilon }$ the symmetrized strain rate tensor. The viscosity

(4)$$\eta = {A\over 2}^{-\lpar {1}/{n}\rpar }\lpar \dot{\epsilon}_{\rm II} + \dot{\epsilon}_0\rpar ^{1-{1}/{n}}$$

is dependent on the second invariant of the strain rate tensor $\dot {\epsilon }_{\rm II}$. Note that we make an isothermal approximation, and take the ice softness parameter $A$ to be a constant. We explicitly note that this assumption may be questionable. However, because models of Greenland thermal conditions frequently do not match borehole observations in the region considered here (e.g. Harrington and others, Reference Harrington, Humphrey and Harper2015) and sliding in this region is an order of magnitude greater than deformation (Maier and others, Reference Maier, Humphrey, Harper and Meierbachtol2019), we choose to avoid the additional computational expense and uncertainty associated with introducing a thermal model. The exponent in Glen's flow law $n = 3$.

Boundary conditions

At the ice surface $\Gamma _{z_{\rm s}}$ and terminal margin $\Gamma _{\rm T}$ (where the ice thickness is assumed to approximate zero), we impose a no-stress boundary condition

(5)$$\tau' \cdot {\bf n} = {\bf 0}\comma\; $$

where ${\bf n}$ is the outward pointing normal vector, and ${\bf 0}$ is the zero vector.

The remaining lateral boundary $\Gamma _{\rm L}$ is synthetic in the sense that there are no natural physical boundary conditions that should be applied there. Here, we adopt the boundary condition of Papanastasiou and others (Reference Papanastasiou, Malamataris and Ellwood1992), who suggest that the boundary term appearing in the weak form of Eqn (1) (the second term in Eqn (B2)) not be replaced by an arbitrary condition (no stress, e.g.), but rather retained and included as an unknown to be determined as part of the solution procedure. Although this does not lead to a unique solution in the strong form of the differential equation, it does lead to one after discretization with the finite element method. The resulting boundary condition for linear Lagrange finite elements specifies that the curvature of both velocity components vanishes at a point near the boundary which for a sufficiently smooth velocity field outside of the domain approximates a stress free boundary at an infinitely distant location. Griffiths (Reference Griffiths1997) refers to this as the ‘no boundary condition’ and show that it is equivalent to solving a reduced order equation in the neighborhood of the boundary, which for the discretization that we describe below reduces to the solution of the shallow ice approximation.

At the basal boundary $\Gamma _{z_{\rm b}}$ we impose the sliding law

(6)$$\tau'\cdot {\bf n} = -\beta^2 N^p \Vert {\bf u}\Vert _2^{q-1} {\bf u}\comma\; $$

with $\beta ^2$ the basal traction coefficient and ${\bf u}$ the ice velocity, and we use $\Vert \cdot \Vert _2$ to denote the standard $L_2$ norm. We note that this sliding law has some theoretical (Fowler, Reference Fowler1987) and empirical (Budd and others, Reference Budd, Keage and Blundy1979; Bindschadler, Reference Bindschadler1983) support, but does not satisfy Iken's bound (Iken, Reference Iken1981). As such there are alternative sliding laws that may be preferable (e.g. Schoof, Reference Schoof2005). However, we defer a detailed comparison of different sliding laws, and condition this study on Eqn (6) being a reasonable (and numerically stable) approximation to the true subglacial process.

The effective pressure $N$ is given by the ice overburden pressure $P_0 = \rho _i g H$ less the water pressure $P_{\rm w}$

$$N = P_0 - P_{\rm w}.$$

The exponents $p$ and $q$ control the non-linear response of basal shear stress to the effective pressure and velocity (respectively). We note several limiting cases of this sliding law: when $p = q = 1$, we recover the linear Budd law (Budd and others, Reference Budd, Keage and Blundy1979). When $p = 0$, we get the pressure-independent Weertman law (Weertman, Reference Weertman1957). In the limit $q\rightarrow \infty$, we recover a perfectly plastic model of basal stress (e.g. Kamb, Reference Kamb1991).

In practice, we use a re-parameterized version of Eqn (6)

(7)$$\tau_{\rm b} = \gamma^2 \hat{N}^p \Vert \hat{{\bf u}}\Vert ^{q-1} \hat{{\bf u}}\comma\; $$

where ${N}/{{\rm Scale}\lpar N\rpar } = \hat {N}$ is the effective pressure non-dimensionalized by the ice overburden averaged over the model domain, and ${{\bf u}}/{{\rm Scale}\lpar {\bf u}\rpar } = \hat {{\bf u}}$ is similar, with the characteristic scale of ${\bf u}$ taken to be 50 m a$^{-1}$. Thus, the resulting relationship between $\gamma ^2$ (which has units of stress) and $\beta ^2$ is

(8)$$\beta^2 = {\gamma^2\over {\rm Scale}\lpar N\rpar ^p {\rm Scale}\lpar {\bf u}\rpar ^q}.$$

This transformation is helpful because the power law terms on the right-hand side of Eqn (6) can vary by several orders of magnitude, thus requiring that $\beta ^2$ does the same in order to maintain a given characteristic surface velocity. The $\gamma ^2$ parameterization circumvents this scale issue. We take $\gamma ^2$, $p$ and $q$ to be unknown but spatially and temporally constant.

Hydrologic model

In order to predict the effective pressure $N$ on which the sliding law depends, we couple the above ice dynamics model to a hydrologic model that simulates the evolution of the subglacial and englacial storage via fluxes of liquid water through an inefficient linked cavity system and an efficient linked channel system. This model closely follows the model GlaDS (Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013), with some alterations in boundary conditions, discretization and opening rate parameterization.

Over a disjoint subdomain $\bar \Omega _i \subset \bar \Omega \comma\; \, \; \bigcup _{i\in {\cal T}} \bar \Omega _i = \bar {\Omega }$, where ${\cal T}$ is the set of triangles in the finite element mesh, the hydraulic potential $\phi = P_{\rm w} + \rho _{\rm w} g z_{\rm b}$ (with $z_{\rm b}$ the bedrock elevation) evolves according to the parabolic equation

(9)$${e_{\rm v}\over \rho_{\rm w} g} {\partial \phi\over \partial t} + \nabla \cdot {\bf q} - {\cal C} + {\cal O} = m\comma\; $$

where $P_{\rm w}$ is the water pressure, $\rho _{\rm w}$ the density of water, ${\bf q}$ the horizontal flux, ${\cal C}$ the rate at which the cavity system closes (pushing water into the englacial system), ${\cal O}$ the rate at which it opens and $m$ is the recharge rate (either from the surface, basal melt or groundwater). The hydraulic potential is related to the effective pressure by

(10)$$N = \rho_{\rm w} g z_{\rm b} + P_0 - \phi.$$

The horizontal flux is given by the Darcy–Weisbach relation

(11)$${\bf q} = -k_{\rm s} h^{\alpha_{\rm s}} \Vert \nabla \phi\Vert _2^{\beta_{\rm s}-2} \nabla \phi\comma\; $$

a non-linear function of the hydraulic potential, characteristic cavity height $h$, bulk conductivity $k_{\rm s}$ and turbulent exponents $\alpha _{\rm s}$ and $\beta _{\rm s}$.

The average subglacial cavity height $h$ evolves according to

(12)$${\partial h\over \partial t} = {\cal O} - {\cal C}.$$

Here, we model the subgrid-scale glacier bed as self-similar, with bedrock asperity heights modeled with a log-normal distribution:

(13)$$\log h_r \sim {\cal N}\lpar \log \bar{h}_r\comma\; \sigma_h^2\rpar \comma\; $$

and a characteristic ratio $r$ of asperity height to spacing. Thus, the opening rate is given by

(14)$${\cal O} = \int_{0}^\infty {\rm Max}\,\left(\Vert {\bf u}\lpar \varsigma = 1\rpar \Vert _2 r \lpar 1 - {h\over h_r}\rpar \comma\; 0\right)\; P\lpar h_r\rpar \; {\rm d}h_r\comma\; $$

where we use $P\lpar{\cdot}\rpar$ to denote the probability density function. For $\sigma _h^2 = 0$, this expression is equivalent to the standard opening rate used in previous studies (e.g. Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013), albeit reparameterized. However, this implies that once the cavity size reaches $h_r$, then the opening rate becomes zero: for a glacier moving increasingly quickly due to a high water pressure, there is no mechanism for subglacial storage capacity to increase. For $\sigma _h^2\gt 0$, our formulation regularizes the opening rate such that there is ‘always a bigger bump’, but with a diminishing effect away from the modal bump size. Here, we make the somewhat arbitrary choice that $\sigma _h^2 = 1$ and take $\bar {h_r}$ to be a tunable parameter.

The cavity closing rate is given by

(15)$${\cal C} = {2\over n^n} A h \vert N\vert ^{n-1} N.$$

Over a domain edge $\partial \Omega _{ij}$ (the edge falling between subdomains $\bar \Omega _i$ and $\bar \Omega _j$), mass conservation implies that

(16)$${\partial S\over \partial t} + {\partial Q\over \partial s} = {{\rm \Xi} - \Pi\over \rho_{\rm w} L} + m_{\rm c}\comma\; $$

with $S$ the size of a channel occurring along that edge, ${\rm \Xi}$ the opening rate due to turbulent dissipation, $\Pi$ the rate of sensible heat changes due to pressure change and $m_{\rm c}$ the exchange of water with adjacent domains. We adopt the constitutive relations given in Werder and others (Reference Werder, Hewitt, Schoof and Flowers2013) for each of these terms. The channel discharge $Q$ is given by another Darcy–Weisbach relation

(17)$$Q = -k_{\rm c} S^\alpha_{\rm c} \bigg\Vert {\partial \phi\over \partial s} \bigg\Vert _2^{\beta-2} {\partial \phi\over \partial s}\comma\; $$

where $k_{\rm c}$ is a bulk conductivity for the efficient channelized system. The channel size evolves according to

(18)$${\partial S\over \partial t} = {{\rm \Xi} - \Pi\over \rho_i L} - {\cal C}_{\rm c}\comma\; $$

with channel closing rate

(19)$${\cal C}_{\rm c} = {2\over n^n} A S \vert N\vert ^{n-1} N.$$

Substitution of Eqn (18) into Eqn (16) leads to an elliptic equation

(20)$${\partial Q\over \partial s} = {{\rm \Xi} - \Pi\over L} \left({1\over \rho_i} - {1\over \rho_{\rm w}}\right) + m_{\rm c}.$$

The exchange term with the surrounding sheet is given by

(21)$$[ {\bf q}\cdot {\bf n}\rsqb _ + + [ {\bf q} \cdot {\bf n}\rsqb _- = m_{\rm c}\comma\; $$

which states that flux into (or out of) a channel is defined implicitly by the flux balance between the two adjacent sheets.

Boundary conditions

We impose a no-flux boundary condition across boundaries $\Gamma _{\rm T} \cup \Gamma _{\rm L}$ in both the sheet and conduit model:

(22)$${\bf q}\cdot {\bf n} = 0$$

(23)$$Q = 0$$

At first glance, this seems to be a strange choice: how then, does water exit the domain? To account for this, we impose the condition that whenever $\phi \gt \phi _{z_{\rm s}}$, where $\phi _{z_{\rm s}} = \rho _{\rm w} g z_{\rm s}$ is the surface potential, any excess water immediately runs off. Because the margins are thin, and the flux across the lateral boundary is zero, the hydraulic head there quickly rises above the level of the ice surface, and the excess water runs off. This heuristic is necessary to avoid the numerically challenging case when potential gradients would imply an influx boundary condition. With a free flux boundary, the model would produce an artificial influx of water from outside the domain in order to keep channels filled, which is particularly problematic in steep topography. Most of the time, the chosen inequality condition has the practical effect of setting the hydraulic potential on the terminus to atmospheric pressure. We note that a better solution would be to devise a model that allows for unfilled conduits along with an explicit modeling of the subaerial hydrologic system. However, we defer that development to later work and condition the results of this study on the heuristic described above.

In addition to this condition, we also enforce the condition that channels do not form at the margins (i.e. $S = 0$ on $\Gamma _{\rm T} \bigcup \Gamma _{\rm L}$). At the terminus, this ensures that there are no channels with unbounded growth perpendicular to the terminus, and also to ensure that lateral boundaries (where $H\gt 0$) are not used as preferential flow paths.

Large ensemble

In order to construct the training data for ${\cal G}$, we must select the values ${\bf m}_i$ over which ${\cal F}$ should be evaluated. Because all values in ${\bf m}$ are positive, yet we do not wish to bias the dataset toward certain regions of the plausible parameter set over others, we choose to draw ${\bf m}$ from a log-uniform distribution with lower and upper bounds ${\bf b}_{\rm L}$ and ${\bf b}_{\rm U}$:

(25)$$\log_{10}\lpar {\bf m}\rpar \sim {\cal U}\lpar {\rm Bound}_{\rm L}\comma\; {\rm Bound}_{\rm U}\rpar .$$

We refer to this distribution as $P_{\rm em}\lpar {\bf m}\rpar$. The specific values of the bounds are given in Table 1, but in general, parameters vary a few orders of magnitude in either direction from values commonly found in the literature. Note that this distribution is not the prior distribution that we will use for Bayesian inference later on. Rather, it is an extremal bound on what we believe viable parameter values to be. However, the support for the distributions is the same, ensuring that the surrogate model is not allowed to extrapolate.

Table 1. Upper and lower bounds for both the log-uniform distribution used to generate surrogate training examples, as well as the log-beta prior distribution

One viable strategy for obtaining training examples would be to simply draw random samples from $P_{\rm em}\lpar {\bf m}\rpar$, and evaluate the high-fidelity model there. However, because we would like to ensure that there is a sample ‘nearby’ all locations in the feasible parameter space, we instead generate the samples using the quasi-random Sobol sequence (Sobol and others, Reference Sobol, Asotsky, Kreinin and Kucherenko2011), which ensures that the parameter space is optimally filled (the sequence is constructed such that the sum of a function evaluated at these samples converges to the associated integral over the domain as quickly as possible). Although the Sobol sequence is defined over the $k$-dimensional unit hypercube, we transform it into a quasi-random sequence in the space of $P_{\rm em}\lpar {\bf m}\rpar$ using the percent point function.

With this distribution of parameters in hand, we evaluate ${\cal F}$ on each sample ${\bf m}_i$. Using 48 cores, this process took $\sim$4 d for 5000 samples. Note that some parameter combinations never converged, in particular cases where $\gamma ^2$ was too low and the resulting velocity fields were many orders of magnitude higher than observed. We discarded those samples and did not use them in subsequent model training.

Surrogate architecture

Dimensionality reduction

We construct the surrogate model ${\cal G}$ in two stages. In the first stage, we perform a PCA ( Shlens, Reference Shlens2014) to extract a limited set of basis functions that can be combined in linear combination such that they explain a maximal fraction of the variability in the ensemble. Specifically, we compute the eigendecomposition

(26)$${\cal S} = V \Lambda V^T\comma\; $$

where $\Lambda$ is a diagonal matrix of eigenvalues and the columns of $V$ are the eigenvectors of the empirical covariance matrix

(27)$${\cal S} = \sum_{i = 1}^m \omega_{d\comma i} \left[\log_{10} {\cal F}\lpar {\bf m}_i\rpar - \log_{10} \bar{{\cal F}} \right]^2\comma\; $$

with $\omega _d$ a vector of weights such that $\sum _{i = 1}^m \omega _{d\comma i} = 1$ (defined later in Eqn (39)) and

(28)$$\log_{10} \bar{{\cal F}} = \sum_{i = 1}^m \omega_{d\comma i} \log_{10} {\cal F}\lpar {\bf m}_i\rpar .$$

We work with log-velocities due to the large variability in the magnitude of fields that are produced by the high-fidelity model.

The columns of $V$ are an optimal basis for describing the variability in the velocities contained in the model ensemble. They represent dominant model modes (Fig. 2) (in the climate literature, these are often called empirical orthogonal functions). We refer to them as ‘eigenglaciers’ in homage to the equivalently defined ‘eigenfaces’ often employed in facial recognition problems (Sirovich and Kirby, Reference Sirovich and Kirby1987). The diagonal entries of $\Lambda$ represent the variance in the data (once again, here these are a large set of model results) explained by each of these eigenglaciers in descending order. As such, we can simplify the representation of the data by assessing the fraction of the variance in the data still unexplained after representing it with $j$ components

(29)$$f\lpar j\rpar = 1 - {\sum_{i = 1}^j \Lambda_{ii}\over \sum_{i = 1}^m \Lambda_{ii}}.$$

Fig. 2. First 12 basis functions in decreasing order of explained variance for one of 50 bootstrap-sampled ensemble members. The color scale is arbitrary.

We find a cutoff threshold $c$ for the number of eigenglaciers to retain by determining $c = \max _j\in \lcub 1\comma\; \, \ldots \comma\; \, m\rcub \colon f\lpar j\rpar \gt s$. We set $s = 10^{-4}$, which is to say that we retain a sufficient number of basis functions such that we can represent 99.99% of the velocity variability in the model ensemble. For the experiments considered here, $c\approx 50$.

Any velocity field that can be produced by the high-fidelity model can be approximately represented as

(30)$${\cal F}\lpar {\bf m}\rpar \approx \sum_{\,j = 1}^c \lambda_{\,j}\lpar {\bf m}\rpar V_j\comma\; $$

where $V_j$ is the $j$-th eigenglacier, and $\lambda _{j}$ is its coefficient. The (row) vector $\lambda \lpar {\bf m}\rpar$ can thus be thought of as a low-dimensional set of ‘knobs’ that control the recovered model output.

Artificial neural network

Unfortunately, we do not a priori know the mapping $\lambda \lpar {\bf m}\rpar$. In the second stage of surrogate creation, we seek to train a function $\lambda \lpar {\bf m}\semicolon\; \theta \rpar$ with trainable parameters $\theta = \lcub W_l\comma\; \, b_l\comma\; \, \alpha _l\comma$ $\beta _l\; \colon \; l = 1\comma\; \, \ldots \comma\; \, L\rcub$ such that the resulting reconstructed velocity field is as close to the high-fidelity model's output as possible, where $L$ is the number of network blocks (see below). For this task, we use a deep but narrow residual neural network. The architecture of this network is shown in Figure 3. Our choice to use a neural network (as opposed to alternative flexible models like Gaussian process regression and polynomial chaos expansion) was motivated primarily by the relatively high dimensionality of our predictions, for which Gaussian processes and polynomial chaos expansions are not well suited due to the difficulty of modeling cross-covariance.

Fig. 3. Architecture of the neural network used as a surrogate model in this study, consisting of four repetitions of linear transformation, layer normalization, dropout and residual connection, followed by projection into the velocity field space through linear combination of basis functions computed via PCA.

As is common for artificial neural networks, we repeatedly apply a four operation block with input $h_{l-1}$ and output $h_l$. As input to the first block we have our parameter vector, so $h_0 = {\bf m}$. In each block, the first operation is a simple linear transformation

(31)$$\hat{{\bf a}}_l = {\bf h}_{l-1} W_l^T + {\bf b}_l\comma\; $$

where $W_l$ and ${\bf b}_l$ are respectively a learnable weight matrix and bias vector for block $l$. To improve the training efficiency of the neural network, the linear transformation is followed by so-called layer normalization (Ba and others, Reference Ba, Kiros and Hinton2016), which $z$-normalizes then rescales the intermediate quantity $\hat {a}_l$

(32)$${\bf a}_l = \alpha_l {\hat{{\bf a}}_l - \mu_l\over \sigma_l} + \beta_l\comma\; $$

where $\mu _l$ and $\sigma _l$ are the mean and standard deviation of $\hat {{\bf a}}_l$, and $\alpha _l$ and $\beta _l$ are learnable layerwise scaling parameters. Next, in order for the artificial neural network to be able to represent non-linear functions, we apply an activation

(33)$$\hat{{\bf z}}_l = {\rm ReLU}\lpar {\bf a}_l\rpar \comma\; $$

where

(34)$${\rm ReLU}\lpar x\rpar = {\rm Max}\lpar x\comma\; 0\rpar $$

is the rectified linear unit (Glorot and others, Reference Glorot, Bordes and Bengio2011). Next we apply dropout (Srivastava and others, Reference Srivastava, Hinton, Krizhevsky, Sutskever and Salakhutdinov2014), which randomly zeros out elements of the activation vector during each epoch of the training phase

(35)$${\bf z}_l = \hat{{\bf z}}_l \odot R\comma\; $$

where $R$ is a vector of Bernoulli distribution random variables with mean $p$. After training is complete and we seek to evaluate the model, we set each element in $R$ to $p$, which implies that the neural network produces deterministic output with each element of the layer weighted by the probability that it was retained during training. Dropout has been shown to effectively reduce overfitting by preventing complex co-adaptation of weights: by never having guaranteed access to a given value during the training phase, the neural network learns to never rely on a single feature in order to make predictions.

Finally, if dimensions allow (which they do for all but the first and last block), the output of the block is produced by adding its input

(36)$${\bf h}_l = {\bf z}_l + {\bf h}_{l-1}\comma\; $$

a so-called residual connection (He and others, Reference He, Zhang, Ren and Sun2016) which provides a ‘shortcut’ for a given block to learn an identity mapping. This mechanism has been shown to facilitate the training of deep neural networks by allowing an unobstructed flow of gradient information from the right end of the neural network (where the data misfit is defined) to any other layer in the network. For each of these intermediate blocks, we utilize $n_h = 64$ nodes.

At the last block as $l = L$, we have that $\lambda \lpar {\bf m}\rpar = {\bf h}_{L-1} W_L^T + {\bf b}_L$. In this study, $L = 5$. $\lambda \lpar {\bf m}\rpar$ is then mapped to a log-velocity field via $V$, as described above. The complete surrogate model is thus defined as

(37)$${\cal G}\lpar {\bf m}\rpar = 10^{\lambda\lpar {\bf m}\rpar V^T}.$$

Surrogate training

To train this model, we minimize the following objective

(38)$$I\lpar \theta\rpar \propto \sum_{i = 1}^m \sum_{\,j = 1}^{n_p} \omega_{d\comma i} A_j \left[\log_{10} {\cal G}\lpar {\bf m}_i\semicolon\; \theta\rpar _j - \log_{10} {\cal F}\lpar {\bf m}_i\rpar _j\right]^2\comma\; $$

where $A_j$ is the fractional area of the $j$-th gridpoint, and $\omega _{d\comma i}\in [ 0\comma\; \, 1] \comma\; \, \sum _{i = 1}^m \omega _{d\comma i} = 1$ is the weight of the $i$-th training example model error. The former term is necessary because our computational mesh resolution is variable, and if were to simply compute the integral as a sum over gridpoints, we would bias the estimator toward regions with high spatial resolution.

The model above is implemented in pytorch, which provides access to objective function gradients via automatic differentiation (Paszke and others, Reference Paszke, Wallach, Larochelle, Beygelzimer, de Alché-Buc, Fox and Garnett2019). We minimize the objective using the ADAM optimizer (Kingma and Ba, Reference Kingma and Ba2014), which is a variant of minibatch stochastic gradient descent. We use a batch size of 64 input–output pairs (i.e. 64 pairs of parameters and their associated high-fidelity model predictions), and begin with a learning rate of $\eta = 10^{-2}$, that is exponentially decayed by one order of magnitude per 1000 epochs (an epoch being one run through all of the training instances). We run the optimization for 4000 epochs.

The results of the surrogate training are shown in Figure 4. We find that for most instances, the surrogate model produces a velocity field in good agreement with the one produced by the high-fidelity model: in most cases the predicted velocities fall within 20% of the high-fidelity model's predictions. Furthermore, in $\gt$99% of instances the nodally defined high-fidelity model predictions fall within three of the ensemble's standard deviations of its own mean (although these residuals are clearly non-normal). The exception to this occurs in instances where the velocity fields are more than three orders of magnitude greater than observations. Since we intend to use the surrogate for inference and such a velocity field implies that the parameters that created it are unlikely to be consistent with observations anyways, this extreme-value misfit will not influence the inference over glacier model parameters.

Fig. 4. Comparison between emulated velocity field (a) and modeled velocity field (b) for three random instances of ${\bf m}$. Note the different velocity scales for each row. These predictions are out of set: the surrogate model was not trained on these examples, and so is not simply memorizing the training data. (c) Difference between high-fidelity and surrogate modeled speeds, normalized by standard deviation of surrogate model ensemble (a $z$-score), with histogram of the same shown by blue line. (d) Difference between high-fidelity and surrogate modeled speeds, normalized by high-fidelity model speeds.

Bayesian bootstrap aggregation

Neural networks are known to be high-variance models, in the sense that while the high-fidelity model may exhibit a monotonic relationship between input parameters and output velocities, the neural network may exhibit high frequency ‘noise’, similar to that exhibited when fitting high-order polynomials to noisy data. This noise is problematic in that it tends to yield local minima that prohibit optimization and sampling procedures from full exploration of the parameter space. In order to reduce this variance, we introduce Bayesian bootstrap aggregation (Breiman, Reference Breiman1996; Clyde and Lee, Reference Clyde and Lee2001) (so-called bagging), in which we train the surrogate described $B$ times, with the sample weights used in Eqn (38) each time randomly drawn from the Dirichlet distribution

(39)$$\omega_{d\comma i} \sim {\rm Dirichlet}\lpar {\bf 1}\rpar \comma\; $$

where ${\bf 1}$ is a vector of ones with length $m$, the number of training instances.

This procedure yields $B$ independent instances of ${\cal G}$ (with single instances hereafter referred to as ${\cal G}_i$), which are combined as a committee. One way to think about this process is that the high-fidelity model is the mean of a distribution, and each ensemble member is a ‘data point’ (a random function) drawn from that distribution. The optimal estimate of the true mean (once again, the high-fidelity model) is the sample mean of the bootstrap samples

(40)$$\bar{{\cal G}}\lpar {\bf m}\rpar = \sum_{i = 1}^B \omega_{e\comma i} {\cal G}_i\lpar {\bf m}\rpar \comma\; $$

with the weights $\omega _{e\comma i}\in [ 0\comma\; \, 1] \comma\; \, \sum _{i = 1}^B\omega _{e\comma i} = 1$. Although this aggregation reduces predictive error (i.e. yields a better approximation to the high-fidelity model) relative to using a single model, uncertainty remains because we are approximating the true mean with the mean based on a finite number of samples. To account for this residual uncertainty in the surrogate model, we can once again employ Bayesian bootstrapping (Rubin, Reference Rubin1981). In principle, we assume that the sample (the computed members of the bagging committee) provide a reasonable approximation to the population (all possible members of the bagging committee) when estimating variability in the mean. In practice, this means that we model $G\lpar {\bf m}\rpar$ as a random function given by Eqn (40) augmented with Dirichlet distributed weights

(41)$$\omega_{e\comma i} \sim {\rm Dirichlet}\lpar {\bf 1}\rpar .$$

Likelihood model

Observations of surface velocity are reported as a field, as are the model predictions, and thus we have an infinite-dimensional Bayesian inference problem (Bui-Thanh and others, Reference Bui-Thanh, Ghattas, Martin and Stadler2013; Petra and others, Reference Petra, Martin, Stadler and Ghattas2014) because there are an infinite number of real-valued coordinates at which to evaluate misfit. However, in contrast to previous studies, rather than finite observations with an infinite parameter space, we have the converse, with continuous (infinite) observations and finite-dimensional parameters. To circumvent this difficulty, we propose a relatively simple approximation that can account for observational correlation and a variable grid size. We first assume a log-likelihood of the form

(44)$$\log P\lpar {\bf d}\vert {\bf m}\comma\; \omega_e\rpar \propto -{1\over 2} \int_{\bar{\Omega}} \int_{\bar{\Omega}'} {r\lpar x\rpar r\lpar x'\rpar \over \sigma\lpar x\comma\; x'\rpar } \rho_d^2 {\rm d}\bar{\Omega}' {\rm d}\bar{\Omega}\comma\; $$

where $\rho _{\rm d}$ is the data density (number of observations per square meter), $\sigma \lpar x\comma\; \, x'\rpar$ is a covariance function

(45)$$\sigma\lpar x\comma\; x'\rpar = \sigma^2_{obs} + \sigma^2_{cor}\left(1 + {d\lpar x\comma\; x'\rpar \over 2 l^2}\right)^{-1}$$

that superimposes white noise with variance $\sigma ^2_{obs}$ and rational exponential noise with variance $\sigma ^2_{cor}$ and characteristic length scale $l$, which we take as four times the local ice thickness. The former term is intended to account for aleatoric observational uncertainty. The latter is a catch-all intended to account for epistemological uncertainty in the flow model and systematic errors in derivation of the velocity fields, with the rational exponential kernel having ‘heavy tails’ that represent our uncertainty in the true correlation length scale of such errors. Although they represent our best efforts at specifying a sensible likelihood model, we emphasize that they are also somewhat arbitrary and can have significant effects on the resulting analysis. However, in the absence of a more clearly justified model, we assume the one presented here.

$r\lpar x\rpar$ is a residual function given by

(46)$$r\lpar x\rpar = \bar{{\cal G}}\lpar x\semicolon\; {\bf m}\comma\; \omega_e\rpar - \Vert {\bf u}_{obs}\Vert _2\lpar x\rpar \comma\; $$

where ${\bf u}_{obs}$ is the satellite derived, annually averaged velocity field described in the ‘Study area’ section, and in which we omit writing the dependence on ${\bf m}$ for brevity.

Because solutions are defined over a finite element mesh, we split the integrals in Eqn (44) into a sum over dual mesh elements $T$ in collection ${\cal T}$

(47)$$\log P\lpar {\bf d}\vert {\bf m}\comma\; \omega_e\rpar \propto -{1\over 2} \sum_{T\in{\cal T}} \sum_{T'\in{\cal T}} \int_T \int_{T'} {r\lpar x\rpar r\lpar x'\rpar \over \sigma\lpar x\comma\; x'\rpar } \rho^2 {\rm d}T' {\rm d}T.$$

Finally, we make the approximation

(48)$$\int_T \int_{T'} {r\lpar x\rpar r\lpar x'\rpar \over \sigma\lpar x\comma\; x'\rpar } \rho^2 {\rm d}T' {\rm d}T \approx {r\lpar x_T\rpar r\lpar x_{T'}\rpar \over \sigma\lpar x_T\comma\; x_{T'}\rpar } \rho^2 A_{T'} A_T\comma\; $$

where $x_T$ are the coordinates of the barycenter of $T$ (the finite element mesh nodes) and $A_T$ its area. Defining

(49)$${\bf r}^T = [ r\lpar x_1\rpar \comma\; r\lpar x_2\rpar \comma\; \ldots\comma\; r\lpar x_N\rpar ] $$

and

(50)$$\Sigma^{-1} = K \hat{\Sigma}^{-1} K\comma\; $$

where $\hat {\Sigma }_{ij} = \sigma \lpar x_i\comma\; \, x_j\rpar$ and $K = {\rm Diag}\lpar [ \rho A_1\comma\; \, \rho A_2\comma\; \, \ldots \comma\; \, \rho A_N] \rpar$ yields the finite-dimensional multivariate-normal likelihood

(51)$$P\lpar {\bf d}\vert {\bf m}\rpar \propto {\rm exp} \left[-{1\over 2} {\bf r}^T \Sigma^{-1} {\bf r} \right].$$

Prior distribution

In principle, we have very little knowledge about the actual values of the parameters that we hope to infer and thus would like to impose a relatively vague prior during the inference process. However, because the surrogate is ignorant of the model physics, we must avoid allowing it to extrapolate beyond the support of the ensemble. One choice that fulfills both of these objectives is to use as a prior the same log-uniform distribution that we used to generate the surrogate. However, the ensemble distribution was designed to cover as broad a support as possible without biasing the surrogate toward fitting parameter values near some kind of mode and does not represent true prior beliefs about the parameter values. Instead, we adopt for the parameters a scaled log-Beta prior

(52)$${\log_{10} {\bf m} - {\rm Bound}_{\rm L}\over {\rm Bound}_{\rm U} - {\rm Bound}_{\rm L}} \sim {\rm Beta}\lpar \alpha = 2\comma\; \beta = 2\rpar $$

This prior reflects our belief that good parameters values are more likely located in the middle of the ensemble, while also ensuring that regions of parameter space outside the support of the ensemble have zero probability.

Marginalization over $\omega _e$

In order to perform the marginalization over bootstrap weights, we make the Monte Carlo approximation

(53)$$\int P\lpar {\bf d}\vert {\bf m}\comma\; \omega_e\rpar P\lpar {\bf m}\rpar P\lpar \omega_e\rpar {\rm d} \omega_e \approx \sum_{i = 1}^N P\lpar {\bf d}\vert {\bf m}\comma\; \omega_{e\comma i}\rpar P\lpar {\bf m}\rpar \comma\; $$

with $\omega _{e\comma i}$ drawn as in Eqn (53), where $N$ is a number of Monte Carlo samples. The terms in the sum are independent, and may be computed in parallel. However, they are also analytically intractable. Thus, we draw samples from each of the summand distributions (the posterior distribution conditioned on an instance of $\omega _e$) using the MCMC procedure described below, then concatenate the sample to form the posterior distribution approximately marginalized over $\omega _e$. The marginalization of the posterior distribution in this way is similar to BayesBag (Bühlmann, Reference Bühlmann2014; Huggins and Miller, Reference Huggins and Miller2019), but with bootstrap sampling applied over models rather than over observations.