{'fraction_non_alphanumeric': 0.03272316444970478, 'fraction_numerical': 0.008137079528183884, 'mean_word_length': 4.913154127564816, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalable and efficient algorithms for the propagation of uncertainty from data through inference to prediction for large-scale problems, with application to flow of the Antarctic ice sheet

The majority of research on efficient and scalable algorithms in computational science and engineering has focused on the forward problem: given parameter inputs, solve the governing equations to determine output quantities of interest. In contrast, here we consider the broader question: given a (large-scale) model containing uncertain parameters, (possibly) noisy observational data, and a prediction quantity of interest, how do we construct efficient and scalable algorithms to (1) infer the model parameters from the data (the deterministic inverse problem), (2) quantify the uncertainty in the inferred parameters (the Bayesian inference problem), and (3) propagate the resulting uncertain parameters through the model to issue predictions with quantified uncertainties (the forward uncertainty propagation problem)? We present efficient and scalable algorithms for this end-to-end, data-to-prediction process under the Gaussian approximation and in the context of modeling the flow of the Antarctic ice sheet and its effect on sea level. The ice is modeled as a viscous, incompressible, creeping, shear-thinning fluid. The observational data come from InSAR satellite measurements of surface ice flow velocity, and the uncertain parameter field to be inferred is the basal sliding parameter. The prediction quantity of interest is the present-day ice mass flux from the Antarctic continent to the ocean. We show that the work required for executing this data-to-prediction process is independent of the state dimension, parameter dimension, data dimension, and number of processor cores. The key to achieving this dimension independence is to exploit the fact that the observational data typically provide only sparse information on model parameters. This property can be exploited to construct a low rank approximation of the linearized parameter-to-observable map.


Introduction
After nearly 2000 years of little change, sea level began to rise in the late 19th century as Earth began to warm, rising 20 cm in the last century. Satellite measurements over the last two decades indicate that this rate is undergoing significant acceleration. While thermal expansion of the oceans and melting of mountain glaciers have been the main contributors, what is of greatest concern is the projected sea level rise from flows of the Greenland and Antarctic ice sheets, as they continue to lose ice to the oceans at increasing rates. Current estimates put this rate at 200 billion metric tons per year; driven by increased warming, collapse of even a small portion of one of these ice sheets has the potential to greatly accelerate this figure. Indeed, recent evidence suggests that sea level rose abruptly at the end of the last interglacial period (118,000 years ago) by 5-6 m; the likely cause is catastrophic collapse of an ice sheet driven by warming oceans [1]. Based on a conservative estimate of a half meter of sea level rise, the Organization for Economic Cooperation and Development estimates that the 136 largest port cities, with 150 million inhabitants and $35 trillion worth of assets, will be at risk from coastal flooding by 2070 [2].
Clearly, model-based projections of the flow of the polar ice sheets will play a central role in anticipating future sea level rise. However, current ice sheet models are subject to considerable uncertainties. Indeed, ice sheet models were left out of the Intergovernmental Panel on Climate Change's 4th Assessment Report, which stated that "the uncertainty in the projections of the land ice contributions [to sea level rise] is dominated by the various uncertainties in the land ice models themselves...rather than in the temperature projections" [3].
Ice is modeled as a three-dimensional, creeping, viscous, incompressible, non-Newtonian, shear-thinning fluid with strain-rate-and temperature-dependent viscosity. Severe mathematical and computational challenges place significant barriers on improving predictability of ice sheet flow models. These include complex and very high-aspect ratio geometry, highly nonlinear and anisotropic rheology, extremely ill-conditioned linear and nonlinear algebraic systems that arise upon discretization as a result of widely-varying viscosity and basal sliding parameter, a broad range of length scales (tens of meters to thousands of kilometers), localization phenomena including fracture, and complex sub-basal hydrological processes.
However, the greatest mathematical and computational challenges lie in quantifying the uncertainties in the predictions of the ice sheet models. These models are characterized by unknown parameter fields describing the basal sliding parameter (resistance to sliding at the base of the ice sheet), basal topography, geothermal heat flux, and rheology. While these parameter fields cannot be directly observed, they can be inferred from satellite observations of ice surface velocities, which leads to a severely ill-posed inverse problem whose solution is extremely challenging. Quantifying the uncertainties that result from inference of these ice sheet parameter fields from noisy data can be accomplished via the framework of Bayesian inference. Upon discretization of the unknown infinite dimensional parameter field, the solution of the Bayesian inference problem takes the form of a very high-dimensional posterior probability density function (pdf) that assigns to any candidate set of parameter fields our belief (expressed as a probability) that a member of this candidate set is the "true" parameter field that gave rise to the observed data. Sampling this posterior pdf to compute, for example, the mean and covariance of the parameters presents tremendous challenges, since not only is it high dimensional, but evaluating the pdf at any point in parameter space requires a forward ice sheet flow simulation-and millions of such evaluations may be required to obtain statistics of interest using state-of-the-art Markov chain Monte Carlo methods. Finally, the ice sheet model parameters and their associated uncertainties can be propagated through the ice sheet flow model to yield predictions of not only the mass flux of ice into the ocean, but also the confidence we have in those predictions. This amounts to solving a system of stochastic PDEs, which again is intractable when the PDEs are complex and highly nonlinear and the parameters are high-dimensional due to discretization of an infinite-dimensional field.
In summary, while one can formulate a data-to-prediction framework to quantify uncertainties from data to inferred model parameters to predictions with an underlying model of non-Newtonian ice sheet flow, attempting to execute this framework for the Antarctic ice sheet (or other large-scale complex models) is intractable for high-dimensional parameter fields using current algorithms. Yet, quantifying the uncertainties in predictions of ice sheet models is essential if these models are to play a significant role in projections of future sea level. The purpose of this paper is to present an integrated framework and efficient, scalable algorithms for carrying out this data-to-prediction process. By scalable, we mean that the cost-measured in number of (linearized) forward (and adjoint) solves-is independent of not only the number of processor cores, but importantly the state variable dimension, the parameter dimension, and the data dimension.
Two key ideas are needed to produce such scalable algorithms. First, we Gaussianize both the posterior pdf that results from Bayesian solution of the inverse problem of inferring ice sheet parameter fields from satellite observations of surface velocity, as well as the pdf resulting from propagating the uncertain parameter fields through the forward ice sheet model to yield predictions of present-day mass flux into the ocean. This is accomplished by linearizing the parameter-to-observable map as well as the parameter-to-prediction map around the maximum a posterior point. We have found that such linearizations for ice sheet flow problems with basal sliding parameter fields to be satisfactory approximations for what would otherwise be an intractable problem [4].
However, even with these linearizations, computing the covariance of each of the resulting pdf's is prohibitive due to the need to solve the forward ice sheet model a number of times equal to the parameter dimension (or data dimension). We overcome this difficulty by recognizing that these maps are inherently low-dimensional, since the data inform a limited number of directions in parameter space, and the predictions are influenced by a limited number of directions in parameter space. Thus, with the right algorithm, the work-as measured by ice sheet Stokes solves-should scale only with the "information dimension." The key idea to achieve this is to construct a low rank approximation of the parameter-to-observable map via a matrix-free randomized SVD method.
The result is a data-to-prediction framework whose computational cost is overwhelmingly dominated by ice sheet model solves, both forward and adjoint (and to a lesser extent elliptic solves representing the action of parameter prior covariances). Scalability of the entire data-to-prediction framework then follows when we show that (1) the number of forward (and adjoint) ice sheet solves needed for the data-to-prediction process is independent of the state dimension, the parameter dimension, and the data dimension; and (2) the forward (and adjoint) ice sheet solver strong-and weak-scales with increasing number of processor cores. We will show that our data-toprediction framework-despite being adamantly "intrusive" to ensure algorithmic scalability of the inversion, uncertainty quantification, and prediction operations-can be expressed in terms of a fixed and dimensionindependent number of forward-like ice sheet model solves, and thus exploits the same algorithms, solver, and parallel implementation needed for the forward problem. Thus, if a scalable (algorithmically and in parallel) forward solver can be designed-as will be shown in Section 2-scalability of the entire data-to-prediction process ensues.
We demonstrate this scalability of the data-to-prediction framework on the problem of predicting the presentday ice mass flux from Antarctica, starting from Interferometric Synthetic Aperture Radar (InSAR) satellite observations of the surface ice flow velocities, inferring the basal sliding parameter field from this data via a 3D nonlinear Stokes ice flow model (Section 3), quantifying the uncertainty in the parameter inference using the Bayesian framework (Section 4), and propagating the uncertain basal sliding parameters through the forward ice sheet flow model to yield predictions of ice mass flux with quantified uncertainties (Section 5).

Forward problem: Modeling ice sheet flow
The forward nonlinear Stokes ice sheet flow solver is the fundamental kernel of our data-to-prediction framework and is invoked repeatedly throughout. It is crucial that this forward solver scales both algorithmically and in parallel. In this section we give a brief summary of the design of our forward ice sheet flow solver and provide performance results. For more details see [5].
The flow of ice is commonly modeled as a non-Newtonian, viscous, incompressible fluid [6,7]. The balance of mass and linear momentum state that where u denotes the ice flow velocity, p the pressure, ρ the mass density of the ice, and g the acceleration of gravity. We employ a constitutive law for ice that relates the stress tensor σ and the strain rate tensoṙ = 1 2 (∇u + ∇u T ) by Glen's flow law [8], where η is the effective viscosity, I is the second order unit tensor,˙ II = 1 2 tr(˙ 2 u ) is the second invariant of the strain rate tensor, n ≥ 1 is Glen's flow law exponent, and A is a flow rate factor that is a function of the temperature T , determined by the Paterson-Budd relation [9]. To construct a temperature field for our simulations, we approximately solve one-dimensional advection-diffusion equations in the vertical direction. The advection velocity for these equations interpolates between the accumulation rate at the upper surface and zero at the base; we use the surface temperature as the upper boundary condition and either the geothermal heat flux or the temperature pressure melting point as the lower boundary condition where appropriate. The surface temperature, accumulation rate, and geothermal heat flux data come from the ALBMAP dataset [10].
The top boundary of the Antarctic ice sheet Γ t is traction-free, and on its bottom boundary Γ b , we impose no-normal flow and a Robin-type condition in the tangential direction: Here, β = β(x ) is the basal sliding parameter field and T := I − n ⊗ n is a projection operator onto the tangential plane, where "⊗" denotes the outer product. Note that β, which relates tangential velocity to tangential traction, subsumes several complex physical phenomena and thus does not itself represent a physical parameter. It depends on the frictional behavior of the ice sheet, the roughness of the bedrock, the thickness of a plastically deforming layer of till, and the amount of water present between the ice sheet and the bedrock. As such it is subject to great uncertainty; subsequent sections will discuss the inverse problem of estimating it from observational data. On lateral boundaries at or below sea-level, we specify that the normal component of stress must match the hydrostatic pressure of sea water; on lateral boundaries above sea-level, we impose no-flow Dirichlet boundary conditions. The lateral boundaries involve neither the observations nor the uncertain basal parameters, and thus in the interests of conciseness, we omit them from the discussion of the deterministic and statistical inverse problems in Sections 3 and 4. The geometric description of the Antarctic ice sheet is constructed from the ALBMAP dataset [10]. We restrict ourselves to the grounded portion of the ice sheet in our simulations, i.e., we neglect ice shelves, the extension of the sheet onto the surface of the ocean. A locally refined mesh of hexahedral elements is used to discretize the ice sheet domain. We construct a coarse quadrilateral mesh that describes the lateral geometry of the ice sheet and use this mesh as the basis for forest-of-quadtree mesh refinement using the p4est library [11]. Each quadrant in the refined mesh is then used as the footprint of a column of hexahedra. We do not constrain the columns to have the same profile, but allow each hexahedron to be independently refined in the vertical direction. This hybrid mesh refinement strategy gives us flexibility to control the quality of the elements in our computational mesh. In particular, we can control the aspect ratio of the elements using local refinement, which is not possible when isotropic refinement is used, such as octree-based refinement.
For accuracy and efficiency, we use a high order accurate and locally volume-conserving discretization that is provably stable for our locally and nonconformingly-refined hexahedral meshes. To be precise, we use the velocity/pressure finite element pair k × disc k−2 for polynomial velocity order k ≥ 2, i.e., with continuous tensorproduct polynomials of order k for each velocity component, and discontinuous tensor-product polynomials of order k − 2 for the pressure [12,13].
An inexact Newton-Krylov method is used to solve the nonlinear systems arising upon discretization of (1), i.e., each Newton linearization is solved inexactly using an iterative Krylov subspace method. If the inexactness is properly controlled, the number of overall Krylov iterations-and thus the overall work-is minimized [5]. Given a velocity/pressure iterate (u j , p j ) at the jth Newton step, the new iterate is computed as (u j+1 , p j+1 ) = (u j , p j ) + α(δ u , δ p ) with an appropriate step length α > 0. The update (δ u , δ p ) T is the solution of the linear system where A(·) corresponds to the linearization of the nonlinear viscous block in (1a) about u j , B is the discretized divergence operator and r j 1 and r j 2 are the nonlinear residuals corresponding to (1a) and (1b), respectively, evaluated at (u j , p j ). We use matrix-free matrix-vector applications to compute linear and nonlinear residuals, which is particularly efficient for high order discretizations as it allows one to exploit the tensor product structure of the finite element basis [14,15,16,5].
Note that due to the form of the nonlinear rheology (1c), the linearization of the nonlinear viscous block in (1a) results in a anisotropic tensor effective viscosity in the operator A(·). Due to self-adjointness, this operator is identical to the viscous block operator of the adjoint Stokes equations. The adjoint equations will be presented in the next section; the expression for the anisotropic tensor effective viscosity is given in (9).
To solve (2) with an iterative Krylov method (GMRES or FGMRES), efficient preconditioning is critical for scalability. We use a block right preconditioner given by where the inverse of the (1,1)-block A(u j ) is approximated by an algebraic multigrid V-cycle (in particular, GAMG from PETSc [17]), andS denotes a Schur complement approximation given by the inverse viscosityweighted lumped mass matrix [18]. Algebraic multigrid for the (1,1)-block requires an assembled fine grid operator, but the matrix A(u j ) is expensive to construct for high order discretizations: for our elements, the assembly time of A(u j ) scales as k 7 . The increasing density of A(u j ) with k also adds to the setup cost of the AMG hierarchy, which is the bottleneck for parallel scalability. We therefore construct an approximationÃ(u j ) for preconditioning based on low order finite elements [16]; in the presence of nonconforming element faces, this can be challenging [5]. The assembly time forÃ(u j ) depends only on the total problem size, not on the order k, and the sparsity is the same as for a low order discretization. Figure 1 defines the Antarctic ice sheet problem we use to study scalability of our nonlinear Stokes solver, described above. The base coarse mesh is shown, from which a sequence of successively finer meshes is constructed. The figure also displays the synthesized basal sliding parameter field β used for the scaling study. The basal sliding parameter field was synthesized using data describing observed surface velocity, thickness, and slope; it is intended to be representative of the "true" β field. Table 1 presents algorithmic and (strong) parallel scaling studies on the sequence of successively refined Antarctic ice flow problems defined in Figure 1. The results demonstrate excellent algorithmic scalability, despite the severe computational challenges of the problem (nonlinearity, indefiniteness, high order discretization, strongly-varying coefficients, anisotropy, high aspect ratio mesh and geometry, locally-refined/nonconforming mesh, multigrid for a vector system). Both the number of Newton iterations and the average number of Krylov iterations per Newton iteration are only mildly sensitive to problem size (over a range from 38M to 2.1B unknowns) and number of cores (from 128 to 131,072 cores). At larger core counts, the setup phase of algebraic multigrid remains the bottleneck in parallel scalability. Geometric multigrid for problems on forest-of-octree meshes has proven efficient and scalable [19], and can eliminate most of the setup time, but must be adapted to the hybrid refinement scheme discussed above. This is a focus of our ongoing work. Now that we have at our disposal a scalable forward nonlinear Stokes solver, we can proceed to tackle the inverse problem and its associated uncertainty, both of which require thousands of forward Stokes solves, as will be seen in the next two sections.

Inverse problem: Inferring the basal sliding parameter from observed ice surface velocity
In this section, we describe the solution of the inverse problem of inferring basal sliding parameters from satellite surface velocity observations with the ice sheet flow model described in the previous section. Several recent studies have focused on inversion for the basal boundary conditions in ice sheet models. Contributions that use adjoint-based gradient information to solve the deterministic inverse problem include [21,22,23,24,25,26,27]. In [28] the authors use InSAR surface velocity measurements, as employed in this section, and invert for the basal sliding parameter field for the Antarctic ice sheet. However, unlike the present work, the authors use the so-called higher-order approximate Stokes model [29,30]. This model is derived from the full Stokes equations by making the assumptions that horizontal gradients of vertical velocities are negligible compared to vertical gradients of horizontal velocities, and that horizontal gradients of vertical shear stresses are small compared to ρ g. The work presented here represents the first continental scale ice sheet inversion using the full Stokes ice flow equations, which is considered the highest fidelity ice flow model available. Moreover, this work is the first to quantify uncertainty in the solution of the Antarctic ice sheet inverse problem using any flow model (Section 4), as well as propagate that uncertainty to a prediction quantity of interest (Section 5).
In the remainder of this section, we present the inverse problem formulation, give a brief overview of the solution method, provide expressions for the gradient based on adjoint ice flow equations, and present results for the deterministic inverse problem.  Table 1: Strong scaling results for our nonlinear Stokes solver for the Antarctic ice sheet problem on ORNL's Titan supercomputer. The table assesses algorithmic and parallel scalability for three problems (P1, P2, and P3) posed on successively finer meshes, each solved to a tolerance of 10 −12 . A description of the problem setup is given in Figure 1. For each run, we report the number of degrees of freedom (#dof), the CPU cores used (#cores), the number of Newton iterations (#Newton), the number of overall preconditioned Krylov iterations (#Krylov), and timings and parallel efficiencies for the matrix-vector multiplications and the multigrid V-cycles (solve time/eff) and for the matrix assembly and the AMG setup (setup time/eff). A block-Jacobi/SOR smoother is used on each level of the V-cycle: incomplete factorization smoothers have been more effective for similar problems [20,5], but were not robust for the (1,1)-block for these problems. For comparison, we show the number of Krylov iterations (#Krylov (Poisson)) necessary to solve a scalar, constant coefficient, linear Poisson problem posed on the same meshes using 3 finite elements and a V-cycle smoothed by damped block-Jacobi/IC(0), which yields closeto-optimal algorithmic scalability. The solution of this highly nonlinear, indefinite, highly heterogeneous, anisotropic, vector-valued ice flow problem requires just ∼ 2× more Krylov iterations per Newton step than is needed for the (linear, positive definite, homogeneous, isotropic, scalar-valued) Poisson problem, indicating just how efficient our algorithms are. The remaining bottleneck for scalability to very large core counts is the setup phase of the algebraic multigrid preconditioner. The " " symbol in the table indicates that the AMG setup times increased with strong scaling, and thus parallel efficiencies are not reported. Figure 1: Left: the coarse quadrilateral mesh from which the 3D meshes used in problems P1, P2, and P3 are extruded and refined. P1 is a discretization of the ice sheet with a maximum element width-to-height aspect ratio of 10:1; P2 and P3 are created by successive bisection in each direction. Each problem uses 3 × disc 1 finite elements, with 330K, 2.6M, and 21M elements, respectively. This amounts to 38M, 270M, and 2.1B unknowns. Right: the exponential of the basal sliding field β used in the Robin-type basal boundary condition for problems P1, P2, and P3. The pre-exponential A −1/n in (1c) is computed from the Paterson-Budd relation [9], with the temperature varying between −50 • C and −10 • C. For Glen's exponent, we use n = 3, and we take ρ = 910kg/m 3 and g = 9.81m/s 2 .

The regularized inverse problem
The inverse problem is formulated as follows: given (possibly noisy) observational data u obs of the ice surface velocity field, we wish to infer the basal sliding parameter field β(x ) (defined as a coefficient in a Robin boundary condition at the base of the ice sheet) that produces a surface velocity field that best fits the observed velocity. This can be formulated as the variational nonlinear least squares optimization problem where the forward velocity u is the solution of the forward nonlinear Stokes problem (1a)-(1d) for a given basal sliding parameter field β, and is an observation operator that restricts the ice velocity field to the top surface. Since the observed surface flow velocities u obs can vary over three or more orders of magnitude, we normalize the data misfit term by |u obs | 2 + , where is a small constant to prevent division by zero.
The regularization term penalizes oscillatory components of the basal sliding parameter field on the basal surface of the ice, thus restricting the solution to smoothly varying fields. Here β 0 is a reference basal sliding parameter, the differential operator is defined as where γ > 0 is the regularization parameter that controls the strength of the imposed smoothness relative to the data misfit, δ > 0, typically small compared to γ, is added to make the regularization operator invertible, ∇ Γ is the tangential gradient, ∆ Γ is the Laplace-Beltrami operator, and ν denotes the outward unit normal on ∂ Γ b . The regularization operator is thus a positive definite elliptic operator of order 4κ. The need for such a regularization stems from the fact that small wavelength components of the basal sliding parameter field cannot be inferred from surface velocity observations due to the smoothing nature of the map from β to u(β) [31].
In the absence of such a term, the inverse problem is ill-posed, that is, its solution is not unique and is highly sensitive to errors in the observations (for general references on regularization in inverse problems, see e.g., [32,33]). When solving the deterministic inverse problem, a value of κ = 1 2 within the regularization operator (4) (which corresponds to classical H 1 -type Tikhonov regularization) suffices for well-posedness. On the other hand, in the Bayesian inference framework (discussed in Section 4), in which the regularization term (4) reflects prior knowledge on the basal sliding parameter field, higher orders of the elliptic regularization operator are required for well-posedness, depending on the spatial dimension of the domain Ω [34]. For example, sufficient values are κ = 1 2 in one dimension and κ = 1 in two and three dimensions.

Inverse problem solver: Inexact Newton-CG
We solve the optimization problem (3) with an inexact matrix-free Newton-CG method. This amounts to (approximately) solving the linear system that arises at each Newton iteration, by the conjugate gradient method, and then updating β by Here, (β) and (β) are the gradient and Hessian, respectively, of (β), both with respect to β, and α k is a step length chosen by a suitable line search method. The next section presents an efficient method for computation of the gradient via a variational adjoint method; in particular, the gradient is given by expression (7), which depends on solutions of forward and adjoint Stokes problems. Hessian computations are made tractable by recalling that the CG solver does not require the Hessian operator by itself; it requires only the action of the Hessian in a given direction. Operator-free (i.e., matrix-free after discretization) computation of this Hessian action (i.e., Hessian-vector product) is presented in Section 4.3, in particular in the expression (21); it involves second order adjoints but has a structure similar to the gradient computation.
Details of the inexact Newton-CG method are provided in [31] and references therein. A summary of the components of the method is as follows: • The Newton system is solved inexactly by early termination of CG iterations via Eisenstat-Walker (to prevent oversolving) and Steihaug (to avoid negative curvature) criteria. • Preconditioning is effected by the inverse of the (elliptic) regularization operator, which is carried out by multigrid V-cycles. • Globalization is by an Armijo backtracking line search.
• Continuation on the regularization parameter is carried out to warm-start the Newton iterations. Parallel and algorithmic scalability follow as a result of the following properties. Because the dominant components of the method can be expressed as solutions or evaluations of Stokes-like systems and inner products, parallel scalability-that is, maintaining high parallel efficiency as the number of cores increases-is assured whenever a scalable solver for the underlying PDEs is available (as demonstrated in Section 2). The remaining ingredient needed to obtain overall scalability is that the method exhibit algorithmic scalability-that is, the number of iterations does not increase with increasing problem size. This is indeed the case: for a wide class of nonlinear inverse problems, the number of outer Newton iterations and of inner CG iterations is independent of the mesh size and hence parameter dimension, as will be demonstrated in Section 3.4. This is a consequence of the use of a Newton solver, of the compactness of the Hessian of the data misfit term in (3), and the preconditioning by the inverse of the regularization operator so that the resulting preconditioned Hessian is a compact perturbation of the identity, for which Krylov subspace methods exhibit mesh-independent iterations [36].

Gradient computation via a variational adjoint method
The Newton-CG method described in the previous section requires computation of the (infinite dimensional) gradient (β), which is the Fréchet derivative of (β) with respect to β. Using the Lagrangian formulation and variational calculus, we can derive an expression for (β) at an arbitrary point β in parameter space [31]; its strong form is given by where the forward velocity/pressure pair (u, p) satisfies the forward Stokes problem (1), and the adjoint velocity/pressure pair (v, q) satisfies the adjoint Stokes problem, Here, Γ accounts for the pointwise scaling by (|u obs | 2 + ) −1 in (3), and the adjoint stress σ v is given in terms of the adjoint strain rate˙ v by where I is the fourth order identity tensor. The adjoint Stokes problem has several notable properties: (1) its only source term is the misfit between observed and predicted surface velocity on the top boundary; (2) while the forward problem is a nonlinear Stokes problem, the adjoint Stokes problem is linear in the adjoint velocity and pressure, and is characterized by a linearized Stokes operator with a 4th-order tensor anisotropic effective viscosity with anisotropic component that acts in the direction of the forward strain rate˙ u ; and (3) the adjoint Stokes operator is the same operator as the Jacobian that arises in Newton's method for solving the forward (nonlinear) Stokes problem (as a consequence of its self-adjointness). Because of this equivalence, we use the same discretization method and mesh, linear solver, and preconditioner for the adjoint problem as was presented for (a Newton step of) the forward problem in Section 2. Use of the same discretization guarantees so-called gradient consistency [31].
In summary, to compute the gradient (7) for a given β field, we first solve the forward nonlinear Stokes problem (1) for the forward velocity/pressure pair (u, p), and then using this forward pair, we solve the adjoint problem (8) for the adjoint velocity/pressure pair (v, q). Finally, both pairs of fields are used to evaluate the gradient expression (7) at any point on the basal surface Γ b or its boundary ∂ Γ b . The action of the Hessian in a given direction could be computed by directional differences of gradients, using the method just described to compute the gradient (necessitating a nonlinear forward Stokes solve for each gradient evaluation). However, it is more efficient to derive expressions for the Hessian action in infinite-dimensional form and compute with them since only a linearized forward solve is required for each Hessian action; presentation of these expressions will be deferred to Section 4.3, where the Hessian plays an important role in characterizing the uncertainty in the solution of the inverse problem.

Inversion results
Here we assess the performance of the inversion algorithm of the previous section, and present continentalscale inversion results for Antarctica. The observational data are the Antarctic ice sheet surface ice velocities obtained using satellite synthetic aperture radar interferometry (InSAR) [37]. We initialize the Newton iteration with a constant initial basal sliding parameter field and iterate until the norm of the gradient of is sufficiently decreased.
First we study the performance of the inexact Newton-CG method described in the previous section as the parameter and data dimensions that characterize the inverse problem grow. Both of these dimensions are tied to the mesh size, since both the basal sliding parameter and observational surface velocity are treated as continuous fields. To make this study tractable, the inversion domain is limited to the Pine Island Glacier region. Performance is characterized in terms of the number of linear or linearized Stokes systems that must be solved, since these  Table 2: Algorithmic scalability of inverse solver for the ice sheet inverse problem posed on the Pine Island Glacier region of Antarctica. Number of Newton and CG iterations and number of linearized Stokes solves for a sequence of increasingly finer meshes (and hence inversion parameters) for the inexact Newton-CG method. The first column (#sdof) shows the number of degrees of freedom for the state variables; the second column (#pdof) shows the number of degrees of freedom for the inversion parameter field; the third column (#N) reports the number of Newton iterations; columns four (#CG) and five (avgCG) report the total and the average (per Newton iteration) number of CG iterations; and the last column (#Stokes) reports the total number of linear(ized) Stokes solves (from forward, adjoint, and incremental forward and adjoint problems). The Newton iterations are terminated when the norm of the gradient is decreased by a factor of 10 5 . The CG iterations are terminated when the norm of the residual of the Newton system drops below a tolerance that is proportional to the norm of the gradient, per the Eisenstat-Walker criterion corresponding to "Choice 2" in [38, Section 2] (with safeguards as described in Section 2.1 of the same). These results illustrate that the cost of solving the inverse problem by the inexact Newton-CG method, measured by the number of Stokes solves, is independent of the parameter dimension. The cost is also independent of the data dimension, since the surface observational velocity field is refined with decreasing mesh size, just as the basal sliding parameter field is refined.
are the core kernels underlying objective function, gradient, and Hessian-vector product computations, and thus overwhelmingly dominate the run time (the cost of the remaining linear algebra associated with the Newton-CG optimization method is negligible relative to that of the Stokes solves). Table 2 presents algorithmic performance for a sequence of increasingly finer meshes, and hence increasing state, parameter, and data dimensions. The regularization operator for these cases is the Laplacian, i.e., in the operator defined by (5), we take k = 1 2 , with δ = 0. We also set the value of the reference basal sliding parameter field β 0 in (4) to 0. An L-curve criterion is employed to determine an optimal value of the regularization parameter. The L-curve criterion requires solution of several inverse problems with different values of γ. Figure 2 depicts the L-curve for the inverse problem defined by the first row of Table 2. Based on this criterion, the regularization parameter for all tests is taken to be 10 −1 . As can be seen in Table 2, the number of Newton and CG iterations required by the inexact Newton-CG inverse solver is insensitive to the parameter and data dimensions, when scaling from 10K to 1.5M inversion parameters (and 96K to 23M states). Thus, the number of Stokes solves does not increase with increasing inverse problem size, leading to a perfectly scalable inverse solver.  Figure 3 depicts the solution of the ice sheet inverse problem, i.e., the reconstruction of the Antarctic basal sliding parameter field β from InSAR observations of surface ice flow velocity and the nonlinear Stokes ice flow model. Note that the basal sliding parameter field varies over nine orders of magnitude. Low (red) and high (blue) values of β represent low and high resistance to basal sliding and correlate with fast and slow ice flows, respectively. As can be seen, weak resistance to basal sliding extends deep into the interior of the continent. Figure 4 attempts to answer the question of how successful is the inversion in creating a model that explains the data well. The top image portrays the observed surface velocity field over the continent, which varies widely from a few meters per year in the interior of the continent (in dark blue), to several kilometers per year in the fast flowing ice streams and outlet glaciers (in bright red). The bottom image shows the surface velocity field that has been reconstructed from solution of the nonlinear Stokes ice flow model using the inferred basal sliding parameter field. By inspection, we can see that the two surface velocity fields agree well, particularly in fast flow regions. The difference between the two images reflects the ill-posed nature of the inverse problem, for which data noise and model inadequacy can amplify errors in poorly-inferred features of the inversion parameter field. Nevertheless, the inversion is successful in inferring a basal sliding parameter field and resulting ice sheet model that is able to fit the data well, particularly in the critical regions of the ice sheet that impact sea level. Despite this success in fitting the data, ultimately we are still dealing with an ill-posed inverse problem in which the "true solution" cannot be found exactly. The question we must address is: what confidence do we have in the inverse solution we have obtained? The deterministic inverse problem described in this section is not equipped to answer this question. Instead, we turn to the framework of Bayesian inference, which provides a systematic means of quantifying uncertainty in the solution of the inverse problem. In the next section we present algorithms for making the Bayesian framework tractable for large-scale ice sheet inverse problems.

Bayesian quantification of parameter uncertainty: estimating the posterior pdf of the basal sliding parameter field
In this section we tackle the problem of quantifying the uncertainty in the solution of the ice sheet inverse problem discussed in the previous section. We adopt the framework of Bayesian inference [39,40], and in particular its extension to infinite-dimensional inverse problems [34]. To keep our discussion compact, we begin by presenting finite-dimensional expressions (i.e., after discretization of the parameter space) for the Bayesian formulation of the inverse problem; we refer the reader to [4] for elaboration of the infinite-dimensional framework and associated discretization issues for the ice sheet inverse problem. Next we discuss a low-rank-based approximation of the posterior covariance (built on a low-rank approximation of the Hessian of the data misfit) that permit scaling to very large parameter dimensions. We then present expressions that show how the lowrank approximation of this Hessian can be computed for the ice sheet inverse problem, and discuss properties of the Hessian that, in combination with the proposed algorithm, permit computation of the uncertainty in the inverse solution in a fixed number of forward/adjoint solves, independent of the parameter dimension. Finally, the methodology is applied to the large-scale Antarctic ice sheet inverse problem.

Bayesian solution of the inverse problem
In the Bayesian formulation, we state the inverse problem as a problem of statistical inference over the space of uncertain parameters, which are to be inferred from data and a physics model. The resulting solution to the statistical inverse problem is a posterior distribution that assigns to any candidate set of parameter fields our belief (expressed as a probability) that a member of this candidate set is the "true" parameter field that gave rise to the observed data. When discretized, this infinite-dimensional inverse problem leads to a large-scale problem of inference over the discrete parameters β ∈ n , corresponding to degrees of freedom in the parameter field mesh.
The posterior probability distribution combines the prior pdf π prior (β) over the parameter space, which encodes any knowledge or assumptions about the parameter space that we may wish to impose before the data are considered, with a likelihood pdf π like (d obs |β), which explicitly represents the probability that a given set of parameters β might give rise to the observed data d obs ∈ m . Bayes' Theorem then states the posterior pdf explicitly as π post (β|d obs ) ∝ π prior (β)π like (d obs |β).
Note that the infinite-dimensional analog of (10) cannot be stated using pdfs but requires Radon-Nikodym derivatives [34]. For many problems, it is reasonable to choose the prior distribution to be Gaussian. If the parameters represent a spatial discretization of a field, the prior covariance operator usually imposes smoothness on the parameters. This is because rough components of the parameter field are typically not observable from the data and must be determined by the prior to result in a well-posed Bayesian inverse problem. Here, we use elliptic PDE operators to construct the prior covariance, which allows us to build on fast, optimal complexity solvers. More precisely, the prior covariance operator is the inverse of the square of the Laplacian-like operator (5), namely prior := −2 = (−γ∆ Γ + δI) −2 , where γ, δ > 0 control the correlation length and the variance of the prior operator. This choice of prior ensures that prior is a trace-class operator, guaranteeing bounded pointwise variance and a well-posed infinite-dimensional Bayesian inverse problem [34,41].
The difference between the observables predicted by the model and the actual observations d obs (velocity observations on the ice sheet's surface) is due to both measurement and model errors, and is represented by the i.i.d. Gaussian random variable "noise" vector e = f (β) − d obs ∈ m , where f (·) is the (generally nonlinear) operator mapping model parameters to output observables. Evaluation of the parameter-to-observable map f (β) requires solution of the nonlinear Stokes ice flow equations (1) followed by extraction of the surface ice flow velocity field, i.e., it is a discretization of u(β) in (3). Restating Bayes' theorem with Gaussian noise and prior, we obtain the statistical solution of the inverse problem, π post (β), as where β prior is the mean of the prior distribution, Γ prior ∈ n×n is the covariance matrix for the prior that arises upon discretization of prior , and Γ noise ∈ m×m is the covariance matrix for the noise, which takes over the role of the scaling matrix Γ from the deterministic formulation presented in section 3.2. As is clear from the expression (11), despite the choice of Gaussian prior and noise probability distributions, the posterior probability distribution need not be Gaussian, due to the nonlinearity of f (β). The non-Gaussianity of the posterior poses challenges for computing statistics of interest for typical large-scale inverse problems, since π post is often a surface in high dimensions, and evaluating each point on this surface requires the solution of the forward model. Numerical quadrature to compute the mean and covariance matrix, for example, is completely out of the question. The method of choice is Markov chain Monte Carlo (MCMC), which judiciously samples the posterior so that sample statistics can be computed. But the use of MCMC for large-scale inverse problems is still prohibitive for expensive forward problems (such as those governed by the nonlinear PDEs of ice sheet flow) and high dimensional parameter spaces (such as the O(10 5 − 10 6 ) parameters characterizing the basal sliding parameter field), since even for modest numbers of parameters, the number of samples required can be in the millions (see for example the discussion in [4,42]). The need to execute millions of forward ice sheet model simulations is simply not feasible, even with today's multi-petaflops systems.
We are thus led to make a quadratic approximation of the negative log of the posterior (11), which results in a Gaussian approximation of the posterior (β MAP , Γ post ). The mean of this posterior distribution, β MAP , is the parameter vector maximizing the posterior (11), and is known as the maximum a posteriori (MAP) point. It can be found by minimizing the negative log of (11), which amounts to solving the optimization problem (3) (i.e., the deterministic inverse problem) with appropriately weighted norms, i.e., The posterior covariance matrix Γ post is then given by the inverse of the Hessian matrix of at β MAP , namely where the Hessian of the misfit is given by Here F is the Jacobian matrix of the parameter-to-observable map evaluated at β MAP , F * is its adjoint, and W uβ , W βu , and W ββ involve second derivatives of the negative log posterior with respect to the parameters and the states (explicit expressions for the ice sheet inverse problem will be given below). The Gaussian approximation will be accurate when the parameter-to-observable map f (β) behaves nearly linearly over the support of the posterior. This will be the case not only when f (β) is weakly nonlinear, but also for directions in parameter space that are poorly informed by the data (in which case f (β) is approximately constant and thus the prior dominates), as well as directions in parameter space that are strongly informed by the data (in which case the posterior variance is small and thus the linearization of f (β) is accurate over the support of the posterior). For the present work, in which we tackle a massive scale Bayesian inverse problem (∼ 10 6 parameters) governed by large-scale ice sheet flow in Antarctica, all known MCMC methods will be prohibitive, so we cannot compare our Hessian-at-the-MAP-based Gaussian approximation with full sampling of the posterior. However, in [4], we compare this Gaussian approximation with several variants of MCMC sampling of the posterior for a model 2D ice flow problem with a moderate number of parameters (∼100), and conclude that the Hessian-based Gaussian approximation can be appropriate, both as a proposal for MCMC, as well as a stand-alone approximation of the posterior.

Low-rank based posterior covariance approximation
As stated above, the Gaussian approximation of the posterior (11), with covariance matrix Γ post (13) that involves the Hessian of the data misfit evaluated at the MAP point (14), is still intractable. The primary difficulty is that the large parameter dimension n prevents any representation of the posterior covariance Γ post as a dense operator. In particular, the Jacobian of the parameter-to-observable map, F , is formally a dense matrix, and requires n forward PDE solves to construct explicitly. This is intractable when n is large and the forward PDEs are expensive, as in our case of Antarctic ice sheet flow. However, a key feature of the operator F is that its action on a (parameter field-like) vector can be formed by solving a (linearized) forward PDE problem; similarly, the action of its adjoint F * on a (observation-like) vector can be formed by solving a (linearized) adjoint PDE. Moreover, for many ill-posed inverse problems, the Hessian matrix of the data misfit term in (3), given by (14), is a discretization of a compact operator, i.e., its eigenvalues collapse to zero. This can be understood intuitively, since only the modes of the parameter field that strongly influence the observations are present in the dominant spectrum of (14). In many ill-posed inverse problems, numerous modes of the parameter field (for example, highly oscillatory ones) will have negligible effect on the observables. The range space thus is effectively finitedimensional even before discretization (and therefore independent of any mesh), and the eigenvalues decay, often rapidly, to zero. Figure 5 illustrates the rapid spectral decay of the data misfit Hessian stemming from our Antarctic ice sheet inverse problem, demonstrating that just 0.3% of the modes of the basal sliding parameter field can be inferred from the data (on the fine mesh). Next we exploit this low-rank structure and the ability to form matrix-free Hessian actions to construct scalable algorithms to approximate the posterior covariance operator.
Rearranging the expression for Γ post in (13) to factor out Γ prior gives This factorization exposes the prior-preconditioned Hessian of the data misfit, Note thatH misfit is symmetric with respect to the inverse prior-weighted inner product. If a square root of Γ prior is available, one can alternatively use a symmetric preconditioning of the data misfit Hessian by the prior covariance operator [43,41]. This results in a prior-preconditioned data misfit Hessian that is symmetric with respect to the Euclidean inner product.
To construct a low-rank approximation ofH misfit , we consider the generalized eigenvalue problem where Λ = diag(λ i ) ∈ n×n contains the generalized eigenvalues and the columns of W ∈ n×n the generalized eigenvectors. Because Γ prior is symmetric positive definite and H misfit is symmetric positive semi-definite, Λ is non-negative and W is orthogonal with respect to Γ −1 prior , so we can scale W so that W T Γ −1 prior W = I. Defining V := Γ −1 prior W , we can rearrange (17) to get When the generalized eigenvalues {λ i } decay rapidly, we can extract a low-rank approximation ofH misfit by retaining only the r largest eigenvalues and corresponding eigenvectors, Here, W r ∈ n×r contains only the r eigenvectors ofH misfit that correspond to the r largest eigenvalues, which are assembled into the diagonal matrix Λ r = diag(λ i ) ∈ r×r , and V r = Γ −1 prior W r . Here, we use randomized SVD algorithms [44,45] to construct the approximate spectral decomposition. Their performance is comparable to Krylov methods (such as Lanczos) and both require only the action of the Hessian on vectors; we show how to do this for the ice sheet Hessian operator in Section 4.3. However, randomized methods have a significant edge over deterministic methods for large-scale problems, since the required Hessian matrix-vector products are independent of each other, providing asynchronicity and fault tolerance.
Once the low-rank approximation (19) has been constructed, we proceed to obtain the posterior covariance matrix. The Sherman-Morrison-Woodbury formula is employed to perform the inverse in (15), where D r := diag(λ i /(λ i + 1)) ∈ r×r . The last term in this expression captures the error due to truncation in terms of the discarded eigenvalues; this provides a criterion for truncating the spectrum, namely that r is chosen such that λ r is small relative to 1. With this low-rank approximation, the final expression for the approximate posterior covariance follows from (15), Note that (20) expresses the posterior uncertainty as the prior uncertainty less any information gained from the data, filtered through the prior (as a consequence of choosing the prior-preconditioned data misfit Hessian as the operator whose spectrum is truncated). The retained eigenvectors of the prior-preconditioned data misfit Hessian are those modes in parameter space that are simultaneously well-informed by the data and assigned high probability by the prior. Figure 6 displays several of these eigenvectors. The low rank update of the prior covariance in (20), which is based on the low rank approximation (19) of the prior-preconditioned Hessian of the data misfit, has recently been shown to the optimal rank-r update with respect to several important loss functions [46].

The action of the Hessian operator
In this section we show that the Hessian-vector products needed in the randomized low-rank algorithms described above can be computed efficiently and matrix-free. The action of the Hessian operator evaluated at a basal sliding parameter field β in a directionβ is given by : Log-linear plot of spectrum of prior-preconditioned data misfit Hessian for two successively finer parameter/state meshes of the inverse ice sheet problem (12). The curves lie on top of each other, indicating mesh independence of the spectrum. This implies that the dominant eigenvalues/eigenvector of the parameter field-and thus the information content of the data, filtered through the priorare independent of the parameter dimension. Moreover, since the surface observations are refined with the state mesh, the invariance of the spectrum to mesh refinement also implies that the dominant eigenvalues are independent of the data dimension (once the information content of the data is resolved). The low rank approximation captures this dominant, data-informed portion of the spectrum. The eigenvalues are truncated at 0.2.
where (û,p) satisfy the incremental forward Stokes equations with the incremental forward stress and (v,q) satisfy the incremental adjoint Stokes equations where the incremental adjoint stress σv is given by and τû = 2η(u)Ψ˙ û , where We see that the incremental forward and incremental adjoint Stokes equations are linearized versions of their forward and adjoint counterparts, having the same (linearized) operator and differing only in the source terms. Figure 6: Eigenvectors of the prior-preconditioned data misfit Hessian corresponding to the 1st, 2nd, 100th, 200th, 500th, and 4000th eigenvalues. Note that different eigenvectors localize in different regions of Antarctica, and that eigenvectors corresponding to smaller eigenvalues are increasingly more oscillatory (and thus inform smaller length scales of the basal sliding parameter field) but are also increasingly less informed by the data.
Thus, we see that computation of Hessian actions on vectors amount to solution of a pair of forward/adjoint (linearized) Stokes equations, similar to the computation of the gradient. Since the gradient and Hessian action require just one and two linearized Stokes solves, respectively, they are significantly cheaper to evaluate than solving the (nonlinear) forward problem, which typically requires an order of magnitude more linearized Stokes solves.

Scalability of Bayesian solution of the inverse problem
We now discuss the overall scalability of our algorithms for Bayesian solution of the inverse problem. First, we characterize the scalability of construction of the low-rank-based approximate posterior covariance matrix in (20). As stated before, the parameter-to-observable map F cannot be constructed explicitly, since it requires n linearized forward Stokes solves. However, as elaborated above, its action on a vector can be computed by a single linearized Stokes solve, regardless of the number of parameters n and observations m. Similarly, the action of F * on a vector can be computed via a linearized adjoint Stokes solve. Moreover, the prior is usually much cheaper to apply than a forward/adjoint solve (here, it is a scalar elliptic solve on the basal boundary). Therefore, the cost of applyingH misfit to a vector-and thus the per iteration cost of the randomized SVD algorithm-is dominated by a pair of linearized forward and adjoint Stokes solves.
The randomized SVD algorithm requires a number of matrix-vector products proportional to the effective rank r of the matrix [45,44]. Thus, the remaining component to establish scalability of the low-rank approximation ofH misfit is independence of r-and therefore the number of matrix-vector products, and hence Stokes solves-from the parameter dimension n. This is the case when H misfit in (14) is a (discretization of a) compact operator, and when preconditioning by Γ prior does not destroy the spectral decay. This situation is typical for many ill-posed inverse problems, in which the prior is either neutral or of smoothing type. Hence, a low-rank approximation ofH misfit can be made that does not scale with parameter or data dimension, instead depending only on the information content of the data, filtered through the prior. This is indeed the case for our ice sheet inverse problem, as demonstrated in Figure 5.
Once the r eigenpairs defining the low rank approximation have been computed, estimates of uncertainty can be computed by interrogating Γ post in (20) at a cost of just r inner products (which are negligible) plus elliptic solves representing the action of the prior Γ prior on a vector (here carried out with an algebraic multigrid solver and therefore algorithmically scalable). For example, samples can be drawn from the Gaussian defined with a covariance Γ post , a row/column of Γ post can be computed, and the action of Γ post in a given direction can be formed, all at cost that is O(r n) in the number of inner products in addition to the O(n) cost of the multigrid solve. Moreover, the posterior variance field, i.e., the diagonal of Γ post , can be found with O(r n) linear algebra plus O(r) multigrid solves.

Samples from the prior and Gassianized posterior
In this section we provide some results that illustrate the quantification of uncertainty in the solution of the Antarctic ice sheet inverse problem. The number of degrees of freedom was 3,785,889 for each of the discretized state-like variables (state, adjoint, incremental state, and incremental adjoint velocities and pressures) and 409,545 for the uncertain basal sliding parameter field. The problem was solved on 1024 processor cores. We used the Laplacian prior operator with κ = 1, γ = 10, and δ = 10 −5 , and took a zero mean of the basal sliding parameter field, i.e., β 0 = 0. To characterize the noise, we used a diagonal noise covariance matrix Γ noise with entries σ i noise = 0.1( d i obs 2 + ) 1/2 for i = 1, ..., m, with = 10 −9 . That is, σ i noise is (at least) 10% of the velocity magnitude at the ith observational data point.
To compute the MAP point via solution of the optimization problem (12), 213 (outer) Newton iterations were necessary to decrease the nonlinear residual by a factor of 10 3 . In each of these outer Newton iterations, the nonlinear forward Stokes problem has to be solved, for which we use an (inner) Newton method. These inner Newton solves are terminated after the residual is decreased by a factor of 10 8 , which takes an average of 5 Newton iterations, each requiring a linearized Stokes solve. In addition to the nonlinear Stokes solve, each (outer) Newton iteration requires computation of the gradient, which costs one linearized Stokes solve, and inexact solution of the (outer) Newton system (6) using CG, which requires a number of Hessian-vector products. The average number of CG iterations per (outer) Newton step was 239. Each Hessian-vector product requires a pair of linearized Stokes solves, one incremental forward and one incremental adjoint. The tolerance for these incremental Stokes solves was set to 10 −6 . Altogether, finding the MAP point required a total of 107,578 linear(ized) Stokes solves.
We approximated the covariance matrix at the MAP point via a low-rank representation employing 5,000 products of the Hessian matrix with random vectors; hence, the cost of this approximation is 10,000 incremental forward/adjoint (linearized) Stokes solves. Thus, the cost (measured in Stokes solves) of quantifying the uncertainty in the MAP solution is about an order of magnitude less than that of finding the MAP point. Figure 7 depicts the MAP point, and Figure 8 shows samples (of the basal sliding parameter field) from the prior (top row) and from the Gaussian approximation of the posterior (bottom row) pdf. The difference between the two sets of samples reflects the information gained from the data in solving the inverse problem. The differences in the basal sliding parameter field across the posterior samples demonstrate that in the fast ice flow regions there is small variability, while in the center and in West Antarctica, we observe larger variability in the inferred parameters, reflecting uncertainty due to insensitivity of surface velocities to the basal sliding in slow velocity and small velocity gradient regions.
This uncertainty in the inference of the basal sliding parameter field, however, is merely an intermediate quantity. What is of ultimate interest is predictions of output quantities of interest, with associated uncertainties, using the ice sheet model with inferred parameters and their associated uncertainties, which have been computed using the methods of this section. This is the subject of the next section.

Prediction with quantified uncertainty: forward propagation of basal sliding parameter uncertainty to mass flux prediction
Once the inverse problem to infer the unknown basal sliding parameter field from observed surface velocities has been solved, and the uncertainty in this inference quantified through a Gaussian approximation of the posterior (made tractable by a low-rank representation of the prior-preconditioned data misfit Hessian), we are ready to propagate the uncertain basal sliding parameter field through the ice flow model to yield a prediction of our quantity of interest with associated uncertainty. Ultimately our interest is in predicting the ice mass flux to the ocean several decades in the future, under various climate change scenarios. However, this requires a model of the ice sheet as an evolving body (a free boundary problem), ideally coupled to an ocean model. In the present work, we have instead chosen to illustrate our data-to-prediction framework with a simpler quantity of interest Q given by the ice mass flux to the ocean using the steady state ice sheet model employed in the previous sections. Below, we describe scalable algorithms for this final step of our data-to-prediction framework.
Formally, this amounts to solving a system of stochastic PDEs given by the nonlinear Stokes forward model with the uncertain basal sliding parameter described by a Gaussian random field. While the low rank approximation of Section 4.2 has resulted in significant dimensionality reduction (from O(10 6 ) to O(10 3 ), as seen in Figure 5), the effective dimension, ∼4000, is still large in absolute terms. Given that this many modes in parameter space are required to quantify the uncertainty in the basal sliding coefficient parameters, and given the expense of solving the large-scale highly-nonlinear forward ice sheet flow problem, the use of Monte Carlo sampling methods would be prohibitive, since millions of forward solves would likely be required to characterize the statistics of the prediction quantity. Similarly, for the problem we target, the use of polynomial chaos methods would be prohibitive due to the curse of dimensionality that afflicts such methods.
Instead, consistent with our Gaussian approximation of the Bayesian solution of the inverse problem, and our desire to scale to very large parameter dimensions, here we linearize the parameter-to-prediction map at the MAP point, resulting in a Gaussian approximation of the prediction pdf, (Q MAP , Γ prediction ). The mean of this prediction pdf, Q MAP , is computed by solving the forward ice sheet flow model (1) using β MAP as the basal sliding parameter, i.e., the MAP point solution of the inverse problem. The covariance operator, Γ prediction , is found by propagating the covariance of the model parameters (which is given by the inverse of the Hessian evaluated at β MAP , i.e., −1 (β MAP )), through the linearized parameter-to-prediction map, i.e., by where (β MAP ) is the Jacobian of the parameter-to-prediction map, evaluated at the MAP point β MAP , and the Hessian at the MAP, (β MAP ), is defined by its action in a direction by (21), which involves solution of the incremental forward (22) and incremental adjoint (23) Stokes problems.
One of the key ideas to enabling scalability of the prediction-under-uncertainty problem is that the Jacobian of the parameter-to-prediction map (β MAP ) can be determined for each prediction quantity Q by computing the gradient of Q with respect to the parameter field β. In our case, we are interested in using the steady state ice where Γ o is an outflow boundary of interest. The gradient of Q with respect to β evaluated at β MAP can then be found as follows. First, solve the forward problem (1) with basal sliding parameter field given by β MAP . Then, solve an adjoint problem defined for the quantity Q, i.e., where the adjoint stress is given by Since we are using the same ice flow model for the prediction problem that was used to infer β, the adjoint problem (26) for Q resembles the adjoint problem for the regularized data misfit functional in (3), but with a different source term given by the variation of Q with respect to the state variables (u, p). This means that the adjoint Stokes solver can be reused for this purpose (which in turn is the same as the linearized forward Stokes solver). Once the forward velocity u and adjoint velocity v are found, (β MAP ), which is the gradient of Q at β MAP , can be found from the gradient expression (7) without the regularization term, i.e., from Of course if a different flow model is used for the prediction phase (such as one governing a dynamically evolving ice sheet), then a new adjoint equation and gradient expression will have to be derived for that model. Now that (β MAP ) has been found, the next step is to form −1 (β MAP ) * (β MAP ), which could be found by solving a linear system using the preconditioned CG method described in Section 3.2. This would require a number of forward/adjoint incremental Stokes solves equal to the number of CG iterations. Instead, since −1 (β MAP ) is available in the compact form given by (20) (based on the low-rank approximation of the priorpreconditioned Hessian of the data misfit), the product −1 (β MAP ) * (β MAP ) can be formed without having to solve incremental Stokes equations, at the cost of just linear algebra (daxpys and inner products). The final step of the inner product of (β MAP ) with −1 (β MAP ) * (β MAP ) is straightforward. The result is the covariance of the ice mass flux Q, given uncertainties in the satellite observations of surface velocity, uncertainties in the basal sliding parameter field as inferred from the satellite data, and the ice sheet flow model and corresponding prediction quantity Q. For a single Q, the prediction density (Q MAP , Γ prediction ) is univariate. The cost to obtain this pdf, once the low rank-based inverse Hessian (20) has been computed in the Bayesian inference phase, is just a forward solve followed by an (as always, linear) adjoint solve, i.e., at little cost beyond the forward solve. Figure 9 provides results of uncertainty propagation from parameters to predictions following the framework outlined above. Recent work on "inference for prediction" has shown that, for linear inverse problems and linear parameter-to-prediction maps, there is a unique direction in parameter space that influences the prediction quantity of interest [47]. Finding this direction involves identifying parameter modes that are both informed by the observational data and also required for estimating the quantity of interest. We adopt the inference-forprediction (IFP) algorithm presented more generally for multiple quantities of interest in [47]. Denoting the Hessian square root by := −1/2 (β MAP ) and using the linearized parameter-to-prediction map (β MAP ), we compute the eigendecomposition ΨΣ 2 Ψ * of * * (β MAP ) (β MAP ) . Since this operator has rank 1, the eigenvector Ψ corresponding to the only nonzero eigenvalue is given by * * (β MAP ), and the corresponding eigenvalue is Σ 2 = * * (β MAP ) 2 . Following [47], the influential direction for prediction (based on linearizations of the parameter-to-observable map and the parameter-to-prediction map) is given by Using the explicit form of Ψ, the influential direction for prediction in the IFP algorithm thus simplifies to = * * (β MAP )Σ −1/2 = Σ −1/2 (β MAP ) −1 * (β MAP ). This direction, or "mode" of the parameter field, is depicted in the bottom row of Figure 9 for ice mass fluxes over three different portions of the outflow boundary (all of Antarctica, just East Antarctica, and just the Totten Glacier). The top row portrays the gradient of Q with respect to β, computed using (27) for each case; these plots show the regions in parameter space to which Q is most sensitive. Note the differences with the bottom row, which captures not only sensitivity of the prediction quantity to β (i.e., the gradient ), but also uncertainty in the β field (i.e., −1 ). The most uncertain regions are not necessarily the most sensitive, and vice versa. Since, as implied by Figure 8, central and West Antarctica exhibit the largest uncertainties, we expect these regions to play an important role in the parameter modes shown in the bottom row of Figure 9, especially where the sensitivity with respect to the basal sliding parameter is not dominant, e.g., in the bottom left figure. On the other hand, in the center and right figures in the bottom row, focusing on East Antarctica, where the posterior samples suggest lower uncertainties, the modes show mixed uncertainty and sensitivity influence.

Conclusions
We have presented a scalable framework for solving end-to-end, data-to-prediction problems, presented in the context of the flow of the Antarctic ice sheet, and motivated by prediction of its contribution to sea level. We begin with observational data and associated uncertainty, then infer model parameters from the data (Section 3), quantify uncertainties in the inference of the parameters (Section 4), and finally propagate those uncertain model parameters to yield prediction quantity of interest with quantified uncertainties (Section 5). We show that the cost of the entire data-to-prediction framework, when measured in forward or adjoint Stokes solves, is a constant independent of the parameter or data dimension. When combined with a scalable forward solver such as that presented in Section 2, this results in a data-to-prediction framework that is independent of the state dimension and number of processor cores as well. Figure 9: The gradient of the prediction quantity Q with respect to the uncertain basal sliding parameter β (top row) and the directions that jointly maximize uncertainty and sensitivity for the prediction (bottom row) following [47]. For better visualization, in the top row we plot ln(| (β MAP )| + 10 −10 ). The quantities of interest are the mass flux of ice into the ocean for the entire boundary (left column), for all of East Antarctica (center column), and for the Totten Glacier region (right column). The mean and standard deviation of the prediction probability distribution for the three ice mass fluxes are: 1170.83 ± 1080.09, 359.60 ± 1.02, and 71.24 ± 0.30 Gt/a, respectively. This scalability is a consequence of three properties of the data-to-prediction process and their exploitation by our algorithms: (1) when inferring the parameter field, the data are informative about only a low-dimensional subspace within the high-dimensional parameter space, and thus a Newton-CG optimization method, preconditioned by the regularization operator, converges in a number of Newton and CG iterations that is independent of the parameter and data dimensions, depending only on the true information contained in the data; (2) when estimating the uncertainty in the inverse solution, the same property dictates that the Hessian of the data misfit admits a low rank representation, and this can be extracted via randomized SVD in a number of matrix-free Hessian-vector products (each of which requires a pair of incremental forward/adjoint Stokes solves) that is also independent of the data and parameter dimensions, and again depends only on the information contained within the data; and (3) when propagating the inferred parameter uncertainties forward through the ice flow model to yield predictions with quantified uncertainties, the prediction pdf can be formed through the action of the inverse Hessian on the gradient of the prediction quantity with respect to the uncertain parameters, which in turn is found through an additional adjoint Stokes solve. Thus, the entire data-to-prediction process is sensitive only to the true information contained within the data, as opposed to the ostensible data or parameter dimensions.
The uncertainty analysis presented here relies on linearizations of the parameter-to-observable and parameterto-prediction maps, leading to Gaussian approximations of the parameter posterior pdf and prediction quantity of interest pdf. Ultimately one would like to relax these approximations and fully explore the resulting non-Gaussian pdfs; how to do this while retaining scalability for such large-scale complex problems remains an open question, and is a subject of ongoing work.