elsa

An isochronal model for ice-sheet layer tracing

elsa — the Englacial Layer Simulation Architecture — advects a stack of isochronal layers through an ice-sheet domain. Its vertical axis is time: each layer is bounded by two isochrones, and mass does not exchange between neighbouring layers. Vertical numerical diffusion is therefore eliminated by construction. elsa is driven by a host ice-sheet model, which supplies the horizontal velocity, the ice thickness, and the surface and basal mass balance; elsa does not modify any of these fields.

The method follows Born (2017), Born and Robinson (2021) and Rieckh et al. (2024). This is a reimplementation of the Bergen ELSA v2.0 rather than a fork. See Design for the departures from the published scheme, and their motivation.

Note

This is a reimplementation from scratch and is under active development. The library, its public API, the NetCDF output, the restart and the full benchmark suite are all in place. The Yelmox coupling is being wired separately.

The interface

A host model needs one module and three calls:

use elsa

type(elsa_class) :: els

call elsa_init(els,"elsa.nml","elsa",time,time_end,xc,yc,zeta_aa,H_ice,"acx_acy")

call elsa_update(els,time,H_ice,ux,uy,smb,bmb)   ! every host timestep

call elsa_end(els)

elsa_update receives an absolute time, computes its own coupling interval, decides internally whether an update is due, and stores its own previous-step ice thickness. The host therefore calls it unconditionally and manages none of elsa’s internal state. Host fields may be supplied in single or double precision; elsa converts them at the boundary. The stagger argument declares where the host’s velocity samples sit — "acx_acy" for staggered faces (as in Yelmo), or "aa" for cell-centred.

What is new relative to v2.0

Correct upwinding

Each cell face is upwinded on its own velocity. v2.0 selected the direction for both x-faces from the sign of one, which produced a negative entry on the matrix diagonal in convergent cells and caused the solve to diverge. Three stability guards existed only to mask this behaviour; all have been removed. See physics.

No linear solver

The advection is now explicit and sub-stepped, so there is no LIS and no system solve. The only dependency of elsa is fesm-utils. v2.0 required LIS, which the user had to build and locate.

Robust interpolation

Any host grid is mapped onto elsa’s grid through one separable, exact stencil, which is built once at initialization. The parameter grid_factor may take any real value ≥ 1, and is no longer restricted to an integer divisor. See interp.

Restart

A run that is stopped and restarted is bit-identical to the run that never stopped. v2.0 could not restart at all. See greenland.

Benchmarks

Every claim above is checked directly. make check runs the Fortran benchmarks, which assert structural properties and exit non-zero on failure; make validate then runs the Julia analysis, which checks the physics against closed-form solutions and produces the figures.

Benchmark	What it establishes
physics	advection is stable, positive and mass-conserving at any CFL
interp	the grid maps are exact
column	isochrones follow the Nye analytic solution
greenland	elsa runs on a real 3D ice sheet, and restarts exactly

Install

elsa is a configme package:

configme install elsa --only     # elsa + fesm-utils