We present a high-level, differentiable, and composable abstraction for the point evaluation of the solution fields of partial differential equation models. The new functionality, embedded in the Firedrake automated finite element system, enables modellers to easily assimilate point data into their models at the point locations, rather than resorting to extrapolation to a computational mesh. We demonstrate the expressiveness and ease with which more mathematically defensible data assimilation can be performed with examples in the fields of groundwater hydrology and glaciology.

In various geoscience disciplines, modellers seek to estimate fields that are challenging to directly observe using measurements of other related fields. These measurements are often sparse, and it is common practice to first extrapolate these measurements to the grid or mesh used for computations. When this estimation procedure is viewed as a deterministic inverse problem, the extrapolation step is undesirable because the choice of extrapolation method introduces an arbitrary algorithmic degree of freedom that can alter the outcomes. When the estimation procedure is instead viewed through the lens of statistical inference, the extrapolation step is undesirable for the additional reason that it obscures the number of statistically independent measurements that are assimilated and thus makes it impossible to apply statistical goodness-of-fit tests or model selection criteria. The introduction of point evaluation into Firedrake, together with its integration into the automatic differentiation features of the system, greatly facilitates the direct assimilation of point data and thus improved methodology for solving both deterministic and statistical inverse problems.

The Pine Island and Thwaites glaciers are the two largest contributors to sea level rise from Antarctica. Here we examine the influence of basal friction and ice shelf basal melt in determining projected losses. We examine both Weertman and Coulomb friction laws with explicit weakening as the ice thins to flotation, which many friction laws include implicitly via the effective pressure. We find relatively small differences with the choice of friction law (Weertman or Coulomb) but find losses to be highly sensitive to the rate at which the basal traction is reduced as the area upstream of the grounding line thins. Consistent with earlier work on Pine Island Glacier, we find sea level contributions from both glaciers to vary linearly with the melt volume averaged over time and space, with little influence from the spatial or temporal distribution of melt. Based on recent estimates of melt from other studies, our simulations suggest that the combined melt-driven and sea level rise contribution from both glaciers may not exceed 10 cm by 2200, although the uncertainty in model parameters allows for larger increases. We do not include other factors, such as ice shelf breakup, that might increase loss, or factors such as increased accumulation and isostatic uplift that may mitigate loss.

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.