Add an ISMIP7 melt-module calibration test group - #982
Conversation
The ismip7_forcing, ismip7_run and forthcoming ismip7_calibration test groups all remap data from an ISMIP7 polar stereographic grid onto a MALI mesh, and all need the same helpers to do it. Those helpers lived inside ismip7_forcing, so using them from another test group would have meant importing across test-group boundaries. Move them to a new landice framework package, compass/landice/ismip7: tests/ismip7_forcing/ice_sheet_params.py -> ismip7/ice_sheet_params.py tests/ismip7_forcing/create_mapfile.py -> ismip7/mapping.py tests/ismip7_forcing/fracture/remap_utils.py -> ismip7/remap.py This also removes real duplication: extrapolate_source had three identical implementations, one in remap_utils and two more as private methods of ProcessThermalForcing and ProcessRunoff. The two copies are deleted and those steps now call the shared function. The fracture forcing files carry a non-CF-compliant integer-year time coordinate, so remap_utils opened them with decode_times=False while the other two copies used the default. Rather than change behavior for any caller, extrapolate_source grows a decode_times argument defaulting to True, and the fracture steps pass False. No behavior change is intended. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Start a new landice test group for calibrating MALI's sub-shelf melt parameterization against the ISMIP7 Antarctic ice-ocean protocol (Reese et al., Sect. 4.2). This commit adds the shared machinery and the first test case: * terms.py -- the four objective-function terms J1-J4, area-weighted and mesh-agnostic. The upstream toolbox assumes a uniform structured grid, converting cell sums with a single scalar cell area, which is wrong on a 4-20 km variable-resolution mesh where cell area varies by a factor of 25. These take an explicit area array instead, so the same code serves both. * quadratic.py -- the Burgard et al. (2022) local quadratic in Python, both as the reference for checking MALI's Fortran and as what the replication drives. * datasets.py -- one registry of which ocean states feed which term, and where the calibration targets live. * objective.py -- assembling the toolbox's arguments and reducing its output to percentiles. * toolbox/ -- the upstream parameter-selection toolbox, vendored verbatim. PROVENANCE.md records the upstream commit and a SHA256, which __init__.py verifies on import so that an edit or a partial update fails loudly rather than quietly changing published numbers. Note that upstream's newest toolbox tag, param-toolbox-v1, predates the last change to the file, so the commit is pinned rather than the tag. * replication/ -- a test case that reproduces the published 8 km calibration through this code path. It needs no MALI run, and it is the check that the vendored toolbox is being driven correctly. The vendored file must stay byte-for-byte identical to upstream, so it is excluded from flake8 and isort. The replication reproduces the published percentiles exactly: K = 4.75e-5 / 8.5e-5 / 1.375e-4, and the Antarctic-mean draft slope comes out at sin(theta) = 0.0051117, matching the protocol's stated 0.005. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The ais test case is the MALI-mesh calibration itself: remap the ISMIP7 masks and the calibration thermal forcing onto the mesh, run a single-timestep MALI melt diagnostic per ocean state and melt form, verify that melt against an independent Python implementation, aggregate it area-weighted, select the parameter, fit dT_b, and report. Two melt forms are calibrated, selected by config: 'ismip7' (the Burgard local quadratic, calibrating K) and 'ismip6' (the non-local quadratic, calibrating gamma0). With a constant salinity the ISMIP6 form is algebraically identical to the Burgard semi-local form, so that is how the semi-local form is calibrated rather than as a third melt module. Only thermal forcing is remapped. The ISMIP7 ocean forcing already processed for the MALI projections carries no salinity, so a projection has no salinity field to read; calibrating against a spatially varying salinity would tune K for physics the projections cannot run. MALI is run with config_ismip7_melt_salinity_source = 'constant' instead. The design decisions that are easy to undo by accident are asserted rather than commented: the ISMIP7/MALI basin numbering is cross-checked against the existing region mask and fails below 80% agreement; the remapped forcing is required to be finite everywhere; the linearity in the melt parameter that licenses one run per ocean state is measured, not assumed; and the melt expression and the vertical interpolation are each checked against an independent implementation. Also add compass's first unit tests, under tests/, run by pytest in CI. These are unit tests: they run in seconds with no input datasets, no MPAS build and no network access. 75 tests cover the area-weighted terms, the melt formulas (including the local/semi-local degeneracy), the vertical interpolation's four code paths, the dataset registry, the config helpers and the vendored toolbox's integrity. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Add user and developer guide pages and API entries, and index them alongside ismip7_forcing and ismip7_run. The pages give the reasoning behind the choices that would otherwise look arbitrary, and that are easy to undo by accident: why only thermal forcing is remapped and salinity is held constant, why there is one MALI run per ocean state rather than one per parameter value, why the two basin numberings are both written under distinct names, why dT_b is fitted after parameter selection rather than before, why the default J4 weighting includes Dotson when the published one does not, and why the minimised objective cannot be compared between melt forms. Also add unit tests for the objective-function assembly, covering the scaling the one-run-per-state design relies on and the weighting rules. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Running the verification against a real MALI ensemble caught two mistakes, both in compass rather than in MALI. The Python reference for the vertical interpolation applied the freezing-point depth correction only where the draft is below the deepest layer centre. MALI applies it in the other downward-extrapolating branch too, where the layer below the draft is beneath the bed. Without it the two disagreed by up to 3.4 K over 11,000 cells. The correction is now applied in both branches, so the check covers all four code paths instead of excluding one. MALI computes melt only where the ice is floating *and* the cell is connected to the open ocean, leaving TFdraft and the melt at zero elsewhere. Compass treated every floating cell as contributing. For an integral that is harmless, since those cells carry zero melt, but it dilutes the area-weighted basin means that J3 is built from. connectedOceanMask is now written to the melt output and both the verification and the aggregation use MALI's own condition, from a single definition in melt_model.read_run. With these fixed, the melt expression agrees with the Python reference to 6.1e-16 and the linearity in the melt parameter is exact. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The linearity check compared the Python reference against itself, which is trivially satisfied and says nothing about MALI. The claim it is meant to support -- that one MALI run per ocean state suffices, rather than one per (state, parameter) pair -- is about MALI. Add two extra runs of the reference ocean state at 0.5 and 2 times the reference melt parameter, and check that each total melt divided by its multiple gives the same number. Only the 'ismip7' form is scaled this way: its parameter is a namelist option, while the ISMIP6 gamma0 is read from an input file. When no scaled runs are set up the check says so rather than reporting a vacuous pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
MALI computes melt on the geometry at the start of the timestep, but the cellMask written to the output is post-step, like thickness and lowerSurface. Over a single step about 5,000 ice-free cells pick up a trace of ice, so the output mask marks them floating even though MALI never computed melt for them. That made the interpolation check compare against values MALI had not calculated -- the last 1.6 K of disagreement -- and it added zero-melt area to the denominator of the area-weighted basin means that J3 is built from. read_run now derives the floating mask from the mesh file with MALI's own test, intersected with connectedOceanMask, so the geometry the mask and the draft come from is the same one the melt was computed on. With this, all three verification checks pass on the real ensemble: the melt expression agrees with the Python reference to 6.1e-16, the vertical interpolation agrees exactly over all four code paths and all 109,412 melting cells, and melt is exactly proportional to the melt parameter. Also fix the replication, where concatenating per-year aggregates collided on a scalar Time coordinate that differs between the observational files. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
TestingAll of the below was run on Chrysalis, with MALI built from Unit tests
|
| percentile | compass | published |
|---|---|---|
| 5th | 4.75000e-05 | 4.75e-05 |
| 50th | 8.50000e-05 | 8.5e-05 |
| 95th | 1.37500e-04 | 1.375e-04 |
The Antarctic-mean ice-draft slope comes out at sin(theta) = 0.0051117, matching the protocol's stated 0.005.
landice/ismip7_calibration/ais
Full ensemble: 56 single-timestep MALI runs (28 ocean states x 2 melt forms) at 128 tasks, plus 2 extra runs for the linearity check. Each run takes about 2 minutes; the ensemble was run as 8 batches of 7 and completed in about 6 minutes wall clock.
Masks. The remapped ISMIP7 basins agree with the mesh's existing regionCellMasks for 98.7% of all cells and 99.1% of floating cells. The one basin that stands out is the Antarctic Peninsula at 82.5% in the theirs-to-ours direction, which is a genuine difference between the ISMIP6-2022 and IMBIE2-v3 outlines rather than an error, and matches what the prototype found.
Verification against an independent Python implementation, on the present-day climatology, over all 109,412 melting cells:
| check | result |
|---|---|
| melt expression, on MALI's own TFdraft | max relative difference 6.05e-16 |
| vertical interpolation, all four code paths | max difference exactly 0.0 K |
| linearity in the melt parameter, measured in MALI | max relative deviation exactly 0 |
The linearity runs give total melt of 892.7163, 1785.4325 and 3570.8650 Gt/yr at 0.5, 1 and 2 times the reference parameter; dividing by the multiple gives 1785.432511 in every case.
Calibration result, 28 ocean states, full weighting, 100,000 samples:
| melt form | 5th | 50th | 95th |
|---|---|---|---|
| ismip7 (K) | 4.75e-05 | 7.75e-05 | 1.15e-04 |
| ismip6 (gamma0, m/yr) | 6750 | 11000 | 19500 |
dT_b fit. All 16 basins fitted, with the largest correction 1.0 K in basin 9 (Eastern Amundsen), comfortably inside the +/-2 K bound the protocol recommends.
Shelf area. MALI's total ice-shelf area is 1510 against an observed 1538 x 10^3 km^2, a ratio of 0.98, with per-basin ratios between 0.53 and 1.05.
Comparison against the prototype
The prototype that produced the science result reported K = 5.00e-5 / 7.75e-5 / 1.15e-4 on the same mesh, computed with a spatially varying salinity. The median and 95th percentile here match exactly and the 5th is one grid step lower, which is the size of shift expected from moving to a constant salinity. The melting-cell count, 109,412, matches the prototype exactly.
The ISMIP6 numbers do not line up as closely: the prototype reported gamma0 = 7453 / 12200 / 21000 against 6750 / 11000 / 19500 here. Superseded. The gamma0 grid stopped at 30,000 m/yr and was truncating a real tail. With the grid widened it gives 7250 / 13000 / 21500 on this mesh, against the prototype's 7453 / 12200 / 21000. See the More Testing comment below; no science review is needed for this after all.
Bugs this testing caught
Three problems were found by running against real data rather than by inspection, all of them in this PR's code rather than in MALI:
- The Python reference applied the freezing-point depth correction in only one of the two branches where MALI extrapolates downward, disagreeing by up to 3.4 K over 11,000 cells.
pyremap'sRemapperinvokesmpirundirectly, which conflicts with the Slurm allocation; the mask remapping now goes through the framework'sbuild_mapping_file, which uses the configuredparallel_executable.- The
cellMaskwritten to the melt output is post-step, and over one timestep about 5,000 ice-free cells pick up a trace of ice. Treating those as floating compared against values MALI never computed, and added zero-melt area to the denominator of the area-weighted basin means J3 is built from. The contributing cells now come from the initial geometry.
The first and third only became visible once the verification ran against a real ensemble, which is the argument for having the verification step at all.
Not tested
The ismip7_forcing consolidation in the first commit is intended to be behaviour-preserving, but I did not run ismip7_forcing before and after to confirm byte-identical output; it needs the CMIP scenario forcing, which is not staged on this machine. The changes are mechanical (moved modules, updated imports, and one decode_times argument added so that each existing caller keeps its previous behaviour), but a reviewer who can run that test group should confirm it.
GrIS was not exercised at all.
Parameter distributions
The distribution of the optimal K over the 100,000 draws, for the ISMIP7 local form — the one the protocol recommends. It is bimodal: a narrow spike at the mode, K = 5.25e-05, and a broad hump centred near 8-9e-05. That structure is a property of the protocol's objective rather than of MALI: the four terms disagree about the best parameter, so draws dominated by one pair of terms land somewhere different from draws dominated by the other. The prototype found the same shape. The whole distribution sits well inside the parameter grid, which runs to 3.0e-04, with no draws at either end, so nothing is being clipped.
The same for the ISMIP6 non-local form, which with a constant salinity is the Burgard semi-local form. Bimodal again, with a sharp spike near gamma0 = 7250 m/yr and a long tail out to about 25,000. Note the small bar hard against the right edge: 104 of the 100,000 draws (0.10%) are pinned at the grid maximum of 30,000 m/yr, so the grid is slightly too narrow for this form. At 0.1% it cannot move the 95th percentile of 19,500, but the upper bound should be raised before these numbers are used for anything, and it is one candidate explanation for why the ISMIP6 result does not reproduce the prototype as closely as the ISMIP7 result reproduces its counterpart.
The two panels are not comparable to each other by the height or position of the objective, only by the parameter each selects — the objective normalises every term by its own median over the parameter ensemble, so its minimised value means nothing across melt forms.
This comment was written by Claude Code.
|
Here's the corresponding branch to target MALI-Dev: See: |
The melt-diagnostic runs needed a graph partition file distributed alongside the mesh, named by a graph_file_prefix config option. A mesh that has no such file could not be used at all, which is the common case for a newly built mesh. Add a make_graph step that builds graph.info from the mesh and partitions it once for the whole ensemble, and drop the config option. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Two problems surfaced when running the calibration on a newly built mesh. build_mapping_file copied the mesh with shutil.copy, which preserves the mode bits, and then wrote lat/lon fields into the copy. A mesh distributed read-only therefore produced a read-only copy and the step died with a permission error. Copy the contents without the mode bits instead, and clear any stale copy first. This affects ismip7_forcing as well, which uses the same framework function. The basin cross-check needs an ISMIP6-era region mask on the same mesh, which a new mesh does not have, so the check that catches the ISMIP7/MALI off-by-one was simply skipped. Add a second check that needs no extra input: Pine Island and Dotson both drain into ISMIP7 basin 9, the Eastern Amundsen, so the shelf mask and the basin field have to agree about that. Using MALI's 1-based numbering where ISMIP7's 0-based is expected puts them in basin 10 and now fails. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Every part of the calibration is a step, so `compass run` from the test case directory does the whole thing and the generated job script runs exactly that; a user config file is the only thing that has to be written by hand. Compass runs steps serially, so the 58 melt diagnostics take a few hours of wall clock. They do not depend on each other, so note that they can be run from their own step directories in separate jobs when that matters. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A draw that chose the largest parameter on the grid wanted a larger one and could not have it. Enough of those and the distribution piles up against the end of the grid, and its upper percentiles say more about where the grid stops than about the objective. On both meshes tried so far, about 0.1% of draws sat on the gamma0 maximum of 30,000 m/yr, which was too few to move the 95th percentile but should not have gone unnoticed. run_optimisation now reports the fraction of draws at the top of the grid and the calibrate step fails above 0.5%. Raise the default gamma0 maximum to 50,000 m/yr so that nothing is clipped. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
More TestingI reran the calibration on the real initial condition,
Verification is unchanged on the new mesh: melt matches the Python reference to 6.0e-16, the vertical interpolation matches exactly over all 108,504 melting cells, and melt is exactly proportional to The Two smaller fixes the new mesh needed: 
Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly. |
1 participant
This adds a
landice/ismip7_calibrationtest group that calibrates MALI's sub-shelf melt parameterization against the ISMIP7 Antarctic ice-ocean protocol (Reese et al., Sect. 4.2), producing the 5th / 50th / 95th percentiles of the melt parameter the ISMIP7 projections need. It has two test cases:replicationreproduces the published 8 km calibration in about a minute with no MALI run, andaisdoes the calibration on a MALI mesh.Four things worth a reviewer's attention:
MALI-Devsubmodule; point[paths] mpas_modelat that build to run theismip7melt form.ismip7_forcinginto a newcompass/landice/ismip7framework package, which also removes three identical copies ofextrapolate_source. Meant to be behaviour-preserving, but it touchesismip7_forcing.tests/, run bypytestin CI.The objective function is vendored verbatim from
ismip7-antarctic-ocean-forcing, which has no installable package;PROVENANCE.mdrecords the upstream commit and a SHA256 that is verified on import and in CI.Posted by Claude Code on @xylar's behalf. The testing, analysis and wording above are AI-authored; please check them accordingly.