Search

Numerical modeling of ice sheet dynamics is a critical tool for projecting future sea-level rise. Among all the processes responsible for the loss of mass of the ice sheets, enhanced ice discharge triggered by the retreat of marine terminating glaciers is one of the key drivers. Numerical models of ice sheet flow are therefore required to include ice front migration in order to reproduce today's mass loss and be able to predict their future. However, the discontinuous nature of calving poses significant numerical challenge for accurately capturing the motion of the ice front. In this study, we explore different stabilization techniques combined with varying reinitialization strategies to enhance the numerical stability and accuracy of solving the level-set function, which tracks the position of the ice front. Through rigorous testing on an idealized domain with semicircular and a straight-line ice front, including scenarios with diverse front velocities, we assess the performance of these techniques. The findings contribute to advancing our ability to model ice sheet dynamics, specifically calving processes, and provide valuable insights into the most effective strategies for simulating and tracking the motion of the ice front. [ABSTRACT FROM AUTHOR]

Ice-sheet flow models capable of accurately projecting their future mass balance constitute tools to improve flood risk assessment and assist sea-level rise mitigation associated with enhanced ice discharge. Some processes that need to be captured, such as grounding line migration, require high spatial resolution (1 km or better). Conventional ice flow models may need significant computational resources because these models mainly execute Central Processing Units (CPUs), which lack massive parallelism capabilities and feature limited peak memory bandwidth. On the other side of the spectrum, Graphics Processing Units (GPUs) are ideally suited for high spatial resolution as the calculations at every grid point can be performed concurrently by thousands of threads or parallel workers. In this study, we combine GPUs with the pseudo-transient (PT) method, an accelerated iterative and matrix-free solving approach, and investigate its performance for finite elements and unstructured meshes applied to two-dimensional (2-D) models of real glaciers at a regional scale. For both Jakobshavn and Pine Island glacier models, the number of nonlinear PT iterations to converge for a given number of vertices N scales in the order of O(N1.2) or better. We compared the performance of PT CUDA C implementation with a standard finite-element CPU-based implementation using the metric: price and power consumption to performance. The single Tesla V100 GPU is 1.5 times the price of the two Intel Xeon Gold 6140 CPU processors. The power consumption of the PT CUDA C implementation was approximately one-seventh of the standard CPU implementation for the test cases chosen in this study. We expect a minimum speed-up of >1.5 to justify the Tesla V100 GPU price to performance. We report the performance (or the speedup) across glacier configurations for degrees of freedom (DoFs) tested to be >1.5 on a Tesla V100. This study is a first step toward leveraging GPU processing power for accurate polar ice discharge predictions. The insights gained from this study will benefit efforts to diminish spatial resolution constraints at increased computing speed. The increased computing speed will allow running ensembles of ice-sheet flow simulations at the continental scale and at high resolution, previously not possible, enabling quantification of model sensitivity to changes in future climate forcings. These findings will be significantly beneficial for process-oriented and sea-level-projection studies over the coming decades. [ABSTRACT FROM AUTHOR]

Time dependent simulations of ice sheets require two equations to be solved: the mass transport equation, derived from the conservation of mass, and the stress balance equation, derived from the conservation of momentum. The mass transport equation controls the advection of ice from the interior of the ice sheet towards its periphery, thereby changing its geometry. Because it is based on a hyperbolic partial differential equation, a stabilization scheme needs to be employed when solved using the finite element method. Several stabilization schemes exist in the finite element method framework, but their respective accuracy and robustness have not yet been systematically assessed for glaciological applications. Here, we compare classical schemes used in the context of the finite element method: (i) Artificial Diffusion, (ii) Streamline Upwinding, (iii) Streamline Upwind Petrov-Galerkin, (iv) Discontinuous Galerkin, and (v) Flux Corrected Transport. We also look at the stress balance equation, which is responsible for computing the ice velocity that `advects' the ice dowstream. To improve the velocity computation accuracy, the ice sheet modeling community employs several sub-element parameterizations of physical processes at the grounding line, the point where the grounded ice starts to float onto the ocean. Here, we introduce a new sub- element parameterization for the driving stress, the force that drives the ice sheet flow. We analyze the response of each stabilization scheme by running transient simulations forced by ice shelf basal melt. The simulations are based on an idealized ice sheet geometry for which there is no influence of bedrock topography. We also perform transient simulations of the Amundsen Sea Sector, West Antarctica, where real bedrock and surface elevations are employed. In both idealized and real ice sheet experiments, stabilization schemes based on artificial diffusion lead systematically to a bias towards more mass loss in comparison to the other schemes, and therefore, should be avoided or employed with a sufficiently high mesh resolution in the vicinity of the grounding line. We also run diagnostic simulations to assess the accuracy of the driving stress parameterization, which in combination with an adequate parameterization for basal stress, provides improved numerical convergence in ice speed computations and more accurate results. [ABSTRACT FROM AUTHOR]

Understanding future impacts of sea-level rise at the local level is paramount to mitigating its effects. In particular, quantifying the range of sea-level rise outcomes in a probabilistic way enables coastal planners to better adapt strategies, depending on cost and timing. For long-term projections, from present-day to the end of the 21st century, frameworks have been developed that provide such probabilistic projections. They rely on sea-level fingerprints where contributions from different processes are sampled at each individual time step and summed up to create probability distributions of sea-level rise for each desired location. While advantageous, this method does not readily allow for including new physics developed in forward models of each component. For example, couplings and feedbacks between ice sheets, ocean circulation, and solid-Earth uplift cannot easily be represented in such frameworks. Indeed, the main impediment to inclusion of more forward model physics in probabilistic sea-level frameworks is the availability of dynamically computed sea-level fingerprints that can be directly linked to local mass changes. Here, we demonstrate such an approach within the Ice-Sheet and Sea-level System Model (ISSM), where we develop a probabilistic framework that can readily be coupled to forward process models such as those for ice sheets, glacial-isostatic adjustment, hydrology and ocean circulation, among others. Through large scale uncertainty quantification, we demonstrate how this approach enables inclusion of incremental improvements in all forward models and provides fidelity to time-correlated processes. The projection system may readily process input and output quantities that are geodetically consistent with space and terrestrial measurement systems. The approach can also account for numerous improvements in our understanding of sea-level processes. [ABSTRACT FROM AUTHOR]

The thermal state of an ice sheet is an important control on its past and future evolution. Some parts of the ice sheets may be polythermal, leading to discontinuous properties at the cold–temperate transition surface (CTS). These discontinuities require a careful treatment in ice sheet models (ISMs). Additionally, the highly anisotropic geometry of the 3D elements in ice sheet modelling poses a problem for stabilization approaches in advection dominated problems. Here, we present extended enthalpy formulations within the finite-element Ice Sheet System Model (ISSM) that show a better performance to earlier implementations. In a first polythermal-slab experiment, we found that the treatment of the discontinuous conductivities at the CTS with a geometric mean produce more accurate results compared to the arithmetic or harmonic mean. This improvement is particularly efficient when applied to coarse vertical resolutions. In a second ice dome experiment, we find that the numerical solution is sensitive to the choice of stabilization parameters in the well-established Streamline Upwind Petrov–Galerkin (SUPG) method. As standard literature values for the SUPG stabilization parameter are not accounting for the highly anisotropic geometry of the 3D elements in ice sheet modelling, we propose a novel Anisotropic SUPG (ASUPG) formulation. This formulation circumvents the problem of high aspect-ratio by treating the horizontal and vertical directions separately in the stabilization coefficients. The ASUPG method provides accurate results for the thermodynamic equation on geometries with very small aspect ratios like ice sheets. [ABSTRACT FROM AUTHOR]

Accurate projections of the evolution of ice sheets in a changing climate require a fine mesh/grid resolution to correctly capture fundamental physical processes, such as the evolution of the grounding line, the region where grounded ice starts to float. The evolution of the grounding line indeed plays a major role in ice sheet dynamics, as it is a fundamental control on marine ice sheet stability. Numerical modeling of grounding line requires significant computational resources since the accuracy of its position depends on grid or mesh resolution. A technique that improves accuracy with reduced computational cost is the adaptive mesh refinement approach, AMR. We present here the implementation of the AMR technique in the finite element Ice Sheet System Model (ISSM) to simulate grounding line dynamics under two different benchmarks, MISMIP3d and MISMIP+. We test different refinement criteria: (a) distance around grounding line, (b) a posteriori error estimator, the Zienkiewicz-Zhu (ZZ) error estimator, and (c) different combinations of (a) and (b). We find that for MISMIP3d setup, refining 5km around the grounding line, both on grounded and floating ice, is sufficient to produce AMR results similar to the ones obtained with uniformly refined meshes. However, for the MISMIP+ setup, we note that there is a minimum distance of 30km around the grounding line required to produce accurate results. We find this AMR mesh-dependency is linked to the complex bedrock topography of MISMIP+. In both benchmarks, the ZZ error estimator presents high values around the grounding line. Particularly for MISMIP+ setup, the estimator also presents high values in the grounded part of the ice sheet, following the complex shape of the bedrock geometry. This estimator helps guide the refinement procedure such that AMR performance is improved. Our results show that computational time with AMR depends on the required accuracy, but in all cases, it is significantly shorter than for uniformly refined meshes. We conclude that AMR without an associated error estimator should be avoided, especially for real glaciers that have a complex bed geometry. [ABSTRACT FROM AUTHOR]

Subglacial hydrology has a significant influence on ice sheet dynamics, yet remains poorly understood. Complex feedbacks play out between the liquid water and the ice, with constantly changing drainage geometry and flow mechanics. A clear tradition has been established in the subglacial hydrology modeling literature of distinguishing between channelized (efficient) and distributed (inefficient) drainage systems or components. Imposing a distinction that changes the governing physics under different flow regimes, however, may not allow for the full array of drainage characteristics to arise. Here, we present a new subglacial hydrology model: SHaKTI (Subglacial Hydrology and Kinetic Transient Interactions). In this model formulation, a single set of governing equations is applied over the entire domain, with a spatially and temporally varying transmissivity that allows for representation of the wide transition between turbulent and laminar flow, and the geometry of each element is allowed to evolve accordingly to form sheet and channel configurations. The model is implemented as a solution in the Ice Sheet System Model (ISSM). We include steady and transient examples to demonstrate features and capabilities of the model, and we are able to reproduce seasonal behavior of the subglacial water pressure that is consistent with observed seasonal velocity behavior in many Greenland outlet glaciers, supporting the notion that subglacial hydrology may be a key influencer in shaping these patterns. [ABSTRACT FROM AUTHOR]

Paleoclimate proxies are being used in conjunction with ice sheet modeling experiments to determine how the Greenland ice sheet responded to past changes, particularly during the last deglaciation. Although these comparisons have been a critical component in our understanding of the Greenland ice sheet sensitivity to past warming, they often rely on modeling experiments that favor minimizing computational expense over increased model physics. Over Paleoclimate timescales, simulating the thermal structure of the ice sheet has large implications on the modeled ice viscosity, which can feedback onto the basal sliding and ice flow. To accurately capture the thermal field, models often require a high number of vertical layers. This is not the case for the stress balance computation, however, where a high vertical resolution is not necessary. Consequently, since stress balance and thermal equations are generally performed on the same mesh, more time is spent on the stress balance computation than is otherwise necessary. For these reasons, running a higher-order ice sheet model (e.g., Blatter-Pattyn) over timescales equivalent to the paleoclimate record has not been possible without incurring a large computational expense. To mitigate this issue, we propose a method that can be implemented within ice sheet models, whereby the vertical interpolation along the z-axis relies on higher-order polynomials, rather than the traditional linear interpolation. This method is tested within the Ice Sheet System Model (ISSM) using quadratic and cubic finite elements for the vertical interpolation on an idealized case and a realistic Greenland configuration. A transient experiment for the ice thickness evolution of a single dome ice sheet demonstrates improved accuracy using the higher-order vertical interpolation compared to models using the linear vertical interpolation, despite having fewer degrees of freedom. This method is also shown to improve a models ability to capture sharp thermal gradients in an ice sheet particularly close to the bed, when compared to models using a linear vertical interpolation. This is corroborated in a thermal steady-state simulation of the Greenland ice sheet using a higher-order model. In general, we find that using a higherorder vertical interpolation decreases the need for a high number of vertical layers, while dramatically reducing model runtime for transient simulations. Results indicate that when using a higher-order vertical interpolation, runtimes for a transient ice sheet relaxation are upwards of 10 to 57 times faster than using a model which has a linear vertical interpolation, and thus requires a higher number of vertical layers to achieve a similar result in simulated ice volume, basal temperature, and ice divide thickness. The findings suggest that this method will allow higher-order models to be used in studies investigating ice sheet behavior over paleoclimate timescales at a fraction of the computational cost than would otherwise be needed for a model using a linear vertical interpolation. [ABSTRACT FROM AUTHOR]

Earth System Models (ESMs) are becoming increasingly complex, requiring extensive knowledge and experience to deploy and use in an efficient manner. They run on high-performance architectures that are significantly different from the everyday environments that scientists use to pre and post-process results (i.e. MATLAB, Python). This results in models that are hard to use for non specialists, and that are increasingly specific in their application. It also makes them relatively inaccessible to the wider science community, not to mention to the general public. Here, we present a new software/model paradigm that attempts to bridge the gap between the science community and the complexity of ESMs, by developing a new JavaScript Application Program Interface (API) for the Ice Sheet System Model (ISSM). The aforementioned API allows Cryosphere Scientists to run ISSM on the client-side of a webpage, within the JavaScript environment. When combined with a Web server running ISSM (using a Python API), it enables the serving of ISSM computations in an easy and straightforward way. The deep integration and similarities between all the APIs in ISSM (MATLAB, Python, and now JavaScript) significantly shortens and simplifies the turnaround of state-of-the-art science runs and their use by the larger community. We demonstrate our approach via a new Virtual Earth System Laboratory (VESL) Web site. [ABSTRACT FROM AUTHOR]

Within the framework of sea-level rise projections, there is a strong need for hindcast validation of the evolution of polar ice sheets in a way that tightly matches observational records (from radar and altimetry observations mainly). However, the computational requirements for making hindcast reconstructions possible are severe and rely mainly on the evaluation of the adjoint state of transient ice-flow models. Here, we look at the computation of adjoints in the context of the NASA/JPL/UCI Ice Sheet System Model, written in C++ and designed for parallel execution with MPI. We present the adaptations required in the way the software is designed and written but also generic adaptations in the tools facilitating the adjoint computations. We concentrate on the use of operator overloading coupled with the AdjoinableMPI library to achieve the adjoint computation of ISSM. We present a comprehensive approach to 1) carry out type changing through ISSM, hence facilitating operator overloading, 2) bind to external solvers such as MUMPS and GSL-LU and 3) handle MPI-based parallelism to scale the capability. We demonstrate the success of the approach by computing sensitivities of hindcast metrics such as the misfit to observed records of surface altimetry on the North-East Greenland Ice Stream, or the misfit to observed records of surface velocities on Upernavik Glacier, Central West Greenland. We also provide metrics for the scalability of the approach, and the expected performance. This approach has the potential of enabling a new generation of hindcast-validated projections that make full use of the wealth of datasets currently being collected, or already collected in Greenland and Antarctica, such as surface altimetry, surface velocities, and/or gravity measurements. [ABSTRACT FROM AUTHOR]

Identifying fast and robust numerical solvers is a critical issue that needs to be addressed in order to improve projections of polar ice sheets evolving in a changing climate. This work evaluates the impact of using sophisticated numerical solvers for transient ice flow simulations using the NASA-JPL/UCI Ice Sheet System Model (ISSM). We identify optimal numerical solvers by testing them on a commonly used ice flow benchmark test, the Ice Sheet Model Intercomparison Project for Higher-Order ice sheet Models (ISMIP-HOM) Experiment F. Three types of analyses are considered: mass transport, horizontal stress balance, and vertical stress balance. A broad suite of solvers is tested, ranging from direct sparse solvers to preconditioned iterative methods. The results of the fastest solvers for each analysis type are ranked based on their scalability across mesh size for each basal sliding conditions specified in Experiment F. We find that the fastest iterative solvers are ~1.5-100 times faster than the default direct solver used in ISSM with speed-ups improving rapidly with increased mesh resolution. We provide a set of recommendations for users in search of efficient solvers to use for transient ice flow simulations, enabling higher-resolution meshes and faster turnaround time. The end result will be improved transient simulations for short-term, highly resolved forward projections (10-100 year time scale) and also improved long-term paleo-reconstructions using higher-order representation of stresses in the ice. This analysis will also enable a new generation of comprehensive uncertainty quantification assessments of forward sea-level rise projections, which rely heavily on ensemble or sampling approaches that are inherently expensive. [ABSTRACT FROM AUTHOR]

Accurately modeling the contribution of Greenland and Antarctica to sea level rise requires to solve partial differential equations at a high spatial resolution. It is important to test the scalability of existing ice sheet models in order to assess whether they are ready to take advantage of new cluster architectures. In this paper, we discuss the overall scaling of the Ice-sheet and Sea-level System Model (ISSM) applied to the Greenland ice sheet. The model setup used as benchmark problem comprises a variety of modules with different levels of complexity and computational demands. The core builds the so-called stress balance module, which uses the higher-order approximation (or Blatter-Pattyn) of the Stokes equations and a mesh of linear prismatic finite elements to compute the ice flow. We develop a detailed user-oriented, yet low-overhead performance instrumentation tailored to the requirements of earth system models and run scaling tests up to 6 144 MPI processes. The results show that the computation of the Greenland model scales overall well up to 3 072 MPI processes, but is eventually slowed down by matrix assembly, the output handling, and lower-dimensional problems that employ lower numbers of unknowns per MPI process. We also discuss improvements of the scaling and identify further improvements needed for climate research. The instrumented version of ISSM, thus, not only identifies potential performance bottlenecks that were not present at lower core counts but also provides the capability to continually monitor the performance of ISSM code basis. This is of long-term significance as the overall performance of ISSM model depends on the subtle interplay between algorithms, their implementation, underlying libraries, compilers, run-time systems and hardware characteristics, all of which are in a constant state of flux. [ABSTRACT FROM AUTHOR]