Solutions to many important partial differential equations satisfy bounds constraints, but approximations computed by finite element or finite difference methods typically fail to respect the same conditions. Chang and Nakshatrala (Comput Methods Appl Mech Eng 320:287–334, 2017) enforce such bounds in finite element methods through the solution of variational inequalities rather than linear variational problems. Here, we provide a theoretical justification for this method, including higher-order discretizations. We prove an abstract best approximation result for the linear variational inequality and estimates showing that bounds-constrained polynomials provide comparable approximation power to standard spaces. For any unconstrained approximation to a function, there exists a constrained approximation which is comparable in the W^1,p norm. In practice, one cannot efficiently represent and manipulate the entire family of bounds-constrained polynomials, but applying bounds constraints to the coefficients of a polynomial in the Bernstein basis guarantees those constraints on the polynomial. Although our theoretical results do not guaruntee high accuracy for this subset of bounds-constrained polynomials, numerical results indicate optimal orders of accuracy for smooth solutions and sharp resolution of features in convection–diffusion problems, all subject to bounds constraints.

The spatial distribution of ocean-induced melting beneath buttressing ice shelves is often cited as an important factor controlling Antarctica’s sea-level contribution. Using numerical simulations, we investigate the relative sensitivity of grounded-ice loss to the spatial distribution and overall volume of ice-shelf melt over two centuries. Contrary to earlier work, we find only minor sensitivity to melt distribution (<6%), with a linear dependence of ice loss on the total melt. Thus, less complex models that need not reproduce the detailed melt distribution may simplify the projection of future sea level. The linear sensitivity suggests a contribution of up to 5.1 cm from Pine Island Glacier over the next two centuries given anticipated levels of ocean warming, provided its ice shelf does not collapse because of other causes.

We introduce a new software package called 'icepack' for modeling the flow of glaciers and ice sheets. The icepack package is built on the finite element modeling library Firedrake, which uses the Unified Form Language (UFL), a domain-specific language embedded into Python for describing weak forms of partial differential equations. The diagnostic models in icepack are formulated through action principles that are specified in UFL. The components of each action functional can be substituted for different forms of the user's choosing, which makes it easy to experiment with the model physics. The action functional itself can be used to define a solver convergence criterion that is independent of the mesh and requires little tuning on the part of the user. The icepack package includes the 2D shallow ice and shallow stream models. We have also defined a 3D hybrid model based on spectral semi-discretization of the Blatter–Pattyn equations. Finally, icepack includes a Gauss–Newton solver for inverse problems that runs substantially faster than the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method often used in the glaciological literature. The overall design philosophy of icepack is to be as usable as possible for a wide a swath of the glaciological community, including both experts and novices in computational science.

Speedup of Pine Island Glacier over the past several decades has made it Antarctica's largest contributor to sea-level rise. The past speedup is largely due to grounding-line retreat in response to ocean-induced thinning that reduced ice-shelf buttressing. While speeds remained fairly steady from 2009 to late 2017, our Copernicus Sentinel 1A/B-derived velocity data show a >12% speedup over the past 3 years, coincident with a 19-km retreat of the ice shelf. We use an ice-flow model to simulate this loss, finding that accelerated calving can explain the recent speedup, independent of the grounding-line, melt-driven processes responsible for past speedups. If the ice shelf’s rapid retreat continues, it could further destabilize the glacier far sooner than would be expected due to surface- or ocean-melting processes.

A number of important questions concern processes at the margins of ice sheets where multiple components of the Earth system, most crucially ice sheets and oceans, interact. Such processes include thermodynamic interaction at the ice–ocean interface, the impact of meltwater on ice shelf cavity circulation, the impact of basal melting of ice shelves on grounded ice dynamics and ocean controls on iceberg calving. These include fundamentally coupled processes in which feedback mechanisms between ice and ocean play an important role. Some of these mechanisms have major implications for humanity, most notably the impact of retreating marine ice sheets on the global sea level. In order to better quantify these mechanisms using computer models, feedbacks need to be incorporated into the modelling system. To achieve this, ocean and ice dynamic models must be coupled, allowing runtime information sharing between components. We have developed a flexible coupling framework based on existing Earth system coupling technologies. The open-source Framework for Ice Sheet–Ocean Coupling (FISOC) provides a modular approach to coupling, facilitating switching between different ice dynamic and ocean components. FISOC allows fully synchronous coupling, in which both ice and ocean run on the same time step, or semi-synchronous coupling in which the ice dynamic model uses a longer time step. Multiple regridding options are available, and there are multiple methods for coupling the sub-ice-shelf cavity geometry. Thermodynamic coupling may also be activated. We present idealized simulations using FISOC with a Stokes flow ice dynamic model coupled to a regional ocean model. We demonstrate the modularity of FISOC by switching between two different regional ocean models and presenting outputs for both. We demonstrate conservation of mass and other verification steps during evolution of an idealized coupled ice–ocean system, both with and without grounding line movement.

A system of subglacial lakes drained on Thwaites Glacier from 2012-2014. To improve coverage for subsequent drainage events, we extended the elevation and icevelocity time series on Thwaites Glacier through austral winter 2019. These new observations document a second drainage cycle in 2017/18 and identified two new lake systems located in the western tributaries of Thwaites and Haynes glaciers. In situ and satellite velocity observations show temporary < 3% speed fluctuations associated with lake drainages. In agreement with previous studies, these observations suggest that active subglacial hydrology has little influence on thinning and retreat of Thwaites Glacier on decadal to centennial timescales.