ACCESS-AIS3 v0.1 - Initial configuration

[lawrenceabird]

This persistent issue is intended to capture discussion and progress of ACCESS-AIS3 v0.1. The initial plan is detailed below, with major configuration steps included as check-marks to track progress. Discussion of ACCESS-AIS3 v0.1 development should take place by adding comments to the thread below.

The ACCESS-AIS3 v0.1 configuration is developed using Jupyter notebooks, but notebook formatting is intended to allow command-line execution as a configuration script in due course.

Configuration steps:

• Mesh
  • Datasets:
    • MEaSUREs v2 Coastline - Buffered 100 km for ocean
  • Refinement:
    • MEaSUREs v2 Velocity field
    • Hmax = 1500 m on fast-flowing grounded regions
    • Hmax = 5000 m on ice-shelves
    • Hmin = 500 m on ice-shelves
• Parameterise
  • Datasets / parameter values:
    • AntarcticaObsISMIP7-v1.1.nc - ISMIP7 MIPKIT contains various datasets/fields used to initialise ISMIP7 models. Where possible, we should use this dataset so that we're consistent with ISMIP7 efforts. This dataset is available here: https://theghub.org/resources/5317
    • Stress balance:
      • Initial inversion target velocity: MEaSUREs v2
      • Initial rigidity: -20 degC
      • Flow Law: Cuffey; n = 3
      • Initial pressure: Overburden
      • Friction law: Schoof
      • Effective Pressure: Output from Ehrenfeucht et al. (2025). Data are available at: https://zenodo.org/records/12738170
      • Bed topography / Surface elevation / Ice thickness / Ice mask: BedMachine v4
    • Thermal balance:
      • Geothermal heat flow: AQ1
      • Surface temp: RACMO 2.4p1
  • Boundary conditions:
    • Stress balance:
      • Grounded ice-front = Dirichlet (?)
      • Floating ice-front = Neumann (?)
    • Thermal balance:
      • Surface temperature spcs from RACMO 2.4p1
• SSA Floating ice rheology inversion (Sensitivity Analysis)
• SSA Floating ice rheology inversion (L-curve)
• SSA Grounded ice friction inversion (Sensitivity Analysis)
• SSA Grounded ice friction inversion (L-curve)
• Extrude
  • Extrude parameterised model to 3D with 10 layers
• HO Floating ice rheology inversion
• HO thermal simulation
• HO Grounded ice friction inversion
• HO Relaxation (100-years)
  • Datasets
    • Mass balance:
      • SMB: RACMO 2.4p1 Mean SMB
      • Basal melt (floating ice): Davison 2023 (for now). Then ISMIP7...?
  • Requested outputs:
    • Regional output for all IMBIE basins included in MEaSUREs v2
    • Grounding line mass flux
    • Grounded area
    • Floating area
    • VAF
    • SMB
    • ...?

[lawrenceabird]

TODO

• Review mesh definition given updated domain extent to include ocean.
• Move Ehrenfeucht et al. (2024) effective pressure dataset into the CCD/ccdtools. See New Dataset Request: Ehrenfeucht et al. (2024) ccdtools#62. Temporarily available here: /g/data/au88/lb9857/access-ais3/data/Antarctica_SubglacialHydrology.nc
• Move to using MIPKIT (rather than individual data files for model initialization). See New Dataset Request: ISMIP7 MIPKIT ccdtools#61
• md.stressbalance.loadingforce doesn't appear to be marshalling correctly during ho_rheology_floating step.
• Need to update marshalling so initial binaries are written to md.cluster.executiondir, rather than $HOME. Large binaries quickly fill $HOME when working with ARE on Gadi. See Add md.cluster.marshalldir to specify where binary files should be marshalled pyISSM#193

---

NOTES

ARE Instance Requirements

This model configuration uses fields from the RACMO2.4p1 dataset. To support this, we require large ARE instance with ~128GB memory. I recommend using a "custom" size (8 CPU; 128GB). Note that smaller memory may be sufficient, but the "actual" memory requirement has not been tested.

Mesh

The initial mesh development was based on an inherited antarctica.exp file. This domain appears to be largely consistent with the MEaSUREs v2 Antarctica Coastline. To support future ocean/ice coupling and the potential need to transfer ocean model information to the offshore-region around Antarctica (to, for example, average offshore temperatures for use in PICO), we should include ocean in the domain. To support this, I have included a process_domain step that generates a new ais_domain.exp file -- this uses the MEaSUREs v2 coastline and adds a 100 km buffer to include ocean area. The new mesh may not yet be refined/optimised for this new domain extent.

Parameterization

Parameterization may need to be revised based on downstream future issues.

pyISSM Development

Relevant pyISSM development for this ACCESS-AIS3 configuration is underway here: ACCESS-NRI/pyISSM#192

[lawrenceabird]

When performing the SSA rheology inversion, we extract floating ice. However, using md.extract() defines Dirichlet BCs around the entire extent of the extracted model. This includes along the ice-front where Neumann are defined by default.

When updating this to re-define Neumann BCs along the ice-front, the inversion does not converge, exceeding the solver residue.

To do this, I am trying:

    print(f"-- Extracting floating ice only (including ice-front)...")
    # Move ice and ocean mask to elements to capture ice-front
    ice_levelset_elements = pyissm.tools.interp.vertex_to_element(md, md.mask.ice_levelset)
    ocean_levelset_elements = pyissm.tools.interp.vertex_to_element(md, md.mask.ocean_levelset)
    mds = md.extract((ice_levelset_elements < 1) & (ocean_levelset_elements < 1))

    print(f"-- Update BCs to ensure Neumann BCs on ice-front...")
    # Update boundary conditions for Neumann BCs on ice-front.
    # NOTE: This is essential to ensure consistent results between
    # inversion and subsequent stressbalance. Since extract() adds Dirichlet
    # BCs around the entire boundary, by default. This is not consistent
    # with the BCs at the ice-front in the main domain, so make them
    # consistent here.
    iceFront = mds.mask.ice_levelset >= 0
    mds.stressbalance.spcvx[iceFront] = np.nan
    mds.stressbalance.spcvy[iceFront] = np.nan
    mds.stressbalance.spcvz[iceFront] = np.nan

I then use this mask for the sensitivity:

    print(f"-- Defining mask to exclude 0 velocity & non-ice nodes from inversion...")
    mask = (mds.inversion.vel_obs > 0) & (mds.mask.ice_levelset < 0)

For now, I have reverted to:

    print('-- Extracting floating ice only...')
    mask = (md.mask.ocean_levelset < 0) & (md.mask.ice_levelset < 0) # Floating ice
    mds = md.extract(mask)

and

    mask = (mds.inversion.vel_obs > 0)

It would be preferred to have consistent BCs at all stages (i.e. set Neumann BCs in the extracted model), but this isn't working at this point. Continuing with the workflow design/setup for now.