We introduce a continuous Galerkin finite element discretization of the non-hydrostatic Boussinesq approximation of the Navier-Stokes equations, suitable for various applications such as coastal ocean dynamics and ice-ocean interactions, among others. In particular, we introduce a consistent modification of the gravity force term which enhances conservation properties for Galerkin methods without strictly enforcing the divergence-free condition. We show that this modification results in a sharp energy estimate, including both kinetic and potential energy. Additionally, we propose a new, symmetric, tensor-based viscosity operator that is especially suitable for modeling turbulence in stratified flow. The viscosity coefficients are constructed using a residual-based shock-capturing method and the method conserves angular momentum and dissipates kinetic energy. We validate our proposed method through numerical tests and use it to model the ocean circulation and basal melting beneath the ice tongue of the Ryder Glacier and the adjacent Sherard Osborn Fjord in two dimensions on a fully unstructured mesh. Our results compare favorably with a standard numerical ocean model, showing better resolved turbulent flow features and reduced artificial diffusion.

The properties of a block preconditioner that has been successfully used in finite element simulations of large scale ice-sheet flow is examined. The type of preconditioner, based on approximating the Schur complement with the mass matrix scaled by the variable viscosity, is well-known in the context of Stokes flow and has previously been analyzed for other types of non-Newtonian fluids. We adapt the theory to hold for the regularized constitutive (power-law) equation for ice and derive eigenvalue bounds of the preconditioned system for both Picard and Newton linearization using inf-sup stable finite elements. The eigenvalue bounds show that viscosity-scaled preconditioning clusters the eigenvalues well with only a weak dependence on the regularization parameter, while the eigenvalue bounds for the traditional non-viscosity-scaled mass-matrix preconditioner are very sensitive to the same regularization parameter. The results are verified numerically in two experiments using a manufactured solution with low regularity and a simulation of glacier flow. The numerical results further show that the computed eigenvalue bounds for the viscosity-scaled preconditioner are nearly independent of the regularization parameter. Experiments are performed using both Taylor-Hood and MINI elements, which are the common choices for inf-sup stable elements in ice-sheet models. Both elements conform well to the theoretical eigenvalue bounds, with MINI elements being more sensitive to the quality of the meshes used in glacier simulations.

This paper concerns a numerical stabilization method for free-surface ice flow called the free-surface stabilization algorithm (FSSA). In the current study, the FSSA is implemented into the numerical ice-flow software Elmer/Ice and tested on synthetic two-dimensional (2D) glaciers, as well as on the real-world glacier of Midtre Lovénbreen, Svalbard. For the synthetic 2D cases it is found that the FSSA method increases the largest stable time-step size at least by a factor of 5 for the case of a gently sloping ice surface (∼ 3°) and by at least a factor of 2 for cases of moderately to steeply inclined surfaces (∼ 6° to 12°) on a fine mesh. Compared with other means of stabilization, the FSSA is the only one in this study that increases largest stable time-step sizes when used alone. Furthermore, the FSSA method increases the overall accuracy for all surface slopes. The largest stable time-step size is found to be smallest for the case of a low sloping surface, despite having overall smaller velocities. For an Arctic-type glacier, Midtre Lovénbreen, the FSSA method doubles the largest stable time-step size; however, the accuracy is in this case slightly lowered in the deeper parts of the glacier, while it increases near edges. The implication is that the non-FSSA method might be more accurate at predicting glacier thinning, while the FSSA method is more suitable for predicting future glacier extent. A possible application of the larger time-step sizes allowed for by the FSSA is for spin-up simulations, where relatively fast-changing climate data can be incorporated on short timescales, while the slow-changing velocity field is updated over larger timescales.

Sea-level projections depend sensitively on the parameterization used for basal slip in glacier flow models. During slip over rock-beds, ice-bed separation increases with slip velocity and basal water pressure. We present a method for using these variables and measured bed topography to estimate the average bed slope in contact with ice, . Three-dimensional numerical modeling of slip over small areas of former beds has shown that changes in with increasing slip velocity and water pressure mimic changes in basal drag. Computed values of can thus provide the form of the slip law that relates drag to velocity and water pressure, avoiding computationally expensive numerical modeling. The method is applied to 618 sections from four former glacier beds. Results generally show an increase in , and hence inferred basal drag, with slip velocity up to a limiting value, consistent with a regularized Coulomb slip law.

Numerical models for predicting future ice mass loss of the Antarctic and Greenland ice sheets require accurately representing their dynamics. Unfortunately, ice-sheet models suffer from a very strict time-step size constraint, which for higher-order models constitutes a severe bottleneck; in each time step a nonlinear and computationally demanding system of equations has to be solved. In this study, stable time-step sizes are increased for a full-Stokes model by implementing a so-called free-surface stabilization algorithm (FSSA). Previously this stabilization has been used successfully in mantle-convection simulations where a similar viscous-flow problem is solved. By numerical investigation it is demonstrated that instabilities on the very thin domains required for ice-sheet modeling behave differently than on the equal-aspect-ratio domains the stabilization has previously been used on. Despite this, and despite the different material properties of ice, it is shown that it is possible to adapt FSSA to work on idealized ice-sheet domains and increase stable time-step sizes by at least one order of magnitude. The FSSA method presented is deemed accurate, efficient and straightforward to implement into existing ice-sheet solvers. 

Glacier slip is usually described using steady-state sliding laws that relate drag, slip velocity and effective pressure, but where subglacial conditions vary rapidly transient effects may influence slip dynamics. Here we use results from a set of laboratory experiments to examine the transient response of glacier slip over a hard bed to velocity perturbations. The drag and cavity evolution from lab experiments are used to parameterize a rate-and-state drag model that is applied to observations of surface velocity and ice-bed separation from the Greenland ice sheet. The drag model successfully predicts observed lags between changes in ice-bed separation and sliding speed. These lags result from the time (or displacement) required for cavities to evolve from one steady-state condition to another. In comparing drag estimates resulting from applying rate-and-state and steady-state slip laws to transient data, we find the peaks in drag are out of phase. This suggests that in locations where subglacial conditions vary on timescales shorter than those needed for cavity adjustment transient slip processes control basal drag.

Ice-sheet responses to climate warming and associated sea-level rise depend sensitively on the form of the slip law that relates drag at the beds of glaciers to their slip velocity and basal water pressure. Process-based models of glacier slip over idealized, hard (rigid) beds with water-filled cavities yield slip laws in which drag decreases with increasing slip velocity or water pressure (rate-weakening drag). We present results of a process-based, three-dimensional model of glacier slip applied to measured bed topographies. We find that consideration of actual glacier beds eliminates or makes insignificant rate-weakening drag, thereby uniting process-based models of slip with some ice-sheet model parameterizations. Computed slip laws have the same form as those indicated by experiments with ice dragged over deformable till, the other common bed condition. Thus, these results may point to a universal slip law that would simplify and improve estimations of glacier discharges to the oceans.

The morphology of glacier beds is a first-order control on their slip speeds and consequent rates of subglacial erosion. As such, constraining the range of bed shapes expected beneath glaciers will improve estimates of glacier slip speeds. To estimate the variability of subglacial bed morphology, we construct 10 high-resolution (10 cm) digital elevation models of proglacial areas near current glacier margins from point clouds produced through a combination of terrestrial laser scanning and photogrammetry techniques. The proglacial areas are located in the Swiss Alps and the Canadian Rockies and consist of predominantly debris-free bedrock of variable lithology (igneous, sedimentary, and metamorphic). We measure eight different spatial parameters intended to describe bed morphologies generated beneath glaciers. Using probability density functions, Bhattacharyya coefficients, principal component analysis, and Bayesian statistical models we investigate the significance of these spatial parameters. We find that the parameters span similar ranges, but the means and standard deviations of the parameter probability density functions are commonly distinct. These results indicate that glacier flow over bedrock may lead to a convergence toward a common bed morphology. However, distinct properties associated with each location prevent morphologies from being uniform.

We investigate the accuracy and robustness of one of the most common methods used in glaciology for finite element discretization of the oe-Stokes equations: linear equal order finite elements with Galerkin least-squares (GLS) stabilization on anisotropic meshes. Furthermore, we compare the results to other stabilized methods. We find that the vertical velocity component is more sensitive to the choice of GLS stabilization parameter than horizontal velocity. Additionally, the accuracy of the vertical velocity component is especially important since errors in this component can cause ice surface instabilities and propagate into future ice volume predictions. If the element cell size is set to the minimum edge length and the stabilization parameter is allowed to vary non-linearly with viscosity, the GLS stabilization parameter found in literature is a good choice on simple domains. However, near ice margins the standard parameter choice may result in significant oscillations in the vertical component of the surface velocity. For these reasons, other stabilization techniques, in particular the interior penalty method, result in better accuracy and are less sensitive to the choice of stabilization parameter. During this work, we also discovered that the manufactured solutions often used to evaluate errors in glaciology are not reliable due to high artificial surface forces at singularities. We perform our numerical experiments in both FEniCS and Elmer/Ice.

To increase our understanding of how ice sheets and glaciers interact with the climate system, numerical models have become an indispensable tool. However, the complexity of these systems and the natural limitation in computational power is reflected in the simplifications of the represented processes and the spatial and temporal resolution of the models. Whether the effect of these limitations is acceptable or not, can be assessed by theoretical considerations and by validating the output of the models against real world data. Equally important is to verify if the numerical implementation and computational method accurately represent the mathematical description of the processes intended to be simulated. This thesis concerns a set of numerical models used in the field of glaciology, how these are applied and how they relate to other study areas in the same field.

The dynamical flow of glaciers, which can be described by a set of non-linear partial differential equations called the Full Stokes equations, is simulated using the finite element method. To reduce the computational cost of the method significantly, it is common to lower the order of the used elements. This results in a loss of stability of the method, but can be remedied by the use of stabilization methods. By numerically studying different stabilization methods and evaluating their suitability, this work contributes to constraining the values of stabilization parameters to be used in ice sheet simulations. Erroneous choices of parameters can lead to oscillations of surface velocities, which affects the long term behavior of the free-surface ice and as a result can have a negative impact on the accuracy of the simulated mass balance of ice sheets.

The amount of basal sliding is an important component that affects the overall dynamics of the ice. A part of this thesis considers different implementations of the basal impenetrability condition that accompanies basal sliding, and shows that methods used in literature can lead to a difference in velocity of 1% to 5% between the considered methods.

The subglacial hydrological system directly influences the glacier's ability to slide and therefore affects the velocity distribution of the ice. The topology and dominant mode of the hydrological system on the ice sheet scale is, however, ill constrained. A third contribution of this thesis is, using the theory of R-channels to implement a simple numerical model of subglacial water flow, to show the sensitivity of subglacial channels to transient processes and that this limits their possible extent. This insight adds to a cross-disciplinary discussion between the different sub-fields of theoretical, field and paleo-glaciology regarding the characteristics of ice sheet subglacial hydrological systems. In the study, we conclude by emphasizing areas of importance where the sub-fields have yet to unify: the spatial extent of channelized subglacial drainage, to what degree specific processes are connected to geomorphic activity and the differences in spatial and temporal scales.

As a whole, the thesis emphasizes the importance of verification of numerical models but also acknowledges the natural limitations of these to represent complex systems. Focusing on keeping numerical ice sheet and glacier models as transparent as possible will benefit end users and facilitate accurate interpretations of the numerical output so it confidently can be used for scientific purposes.