Create 3d Greenland thermal forcing from ISMIP7-provided 2d and EN4 data - #979
Merged
trhille merged 29 commits intoSep 28, 2026
Merged
Conversation
trhille
force-pushed
the
landice/create_3d_ismip7_greenland_tf
branch
3 times, most recently
from
September 21, 2026 17:19
95b3207 to
a263d2e
Compare
Port the standalone greenland_thermal_forcing tool into a new optional step of the ocean_thermal test case. build_3d_thermal_forcing (GrIS only, gated by process_ocean_thermal_3d) converts the 2-D GrIS thermal forcing into a 30-level 3-D field for MALI's nonlocal (Jourdain et al. 2020) melt scheme: seven regional EN4 vertical profiles, per-cell seafloor anchoring to the 2-D forcing, and per-region deltaT calibration. It auto-chains from the 2-D output the same run just produced and writes ismip6shelfMelt_3dThermalForcing (plus deltaT/gamma0/zOcean/basin), supplementing the 2-D file. The ported science lives in ocean_thermal/greenland_3d.py; the 3-D-specific parameters are supplied via a JSON config (config_file), while mesh, 2-D forcing, output, and diagnostics paths are injected from the compass config. Also standardize ocean scenario output names to 2dThermalForcing (GrIS) / 3dThermalForcing (AIS) and update the ismip7_run tf globs accordingly (AIS ingests 3-D, GrIS 2-D). Adds docs and no new dependencies (shapely/h5py/dask already present). Not yet wired: ismip7_run GrIS streams to ingest the 3-D field (follow-up).
Add use_3d_thermal_forcing (default false) to [ismip7_run_gris]. When true, ismip7_gris looks for *3dThermalForcing_*.nc instead of *2dThermalForcing_*.nc, reads ismip6shelfMelt_3dThermalForcing + ismip6shelfMelt_zOcean at annual intervals (instead of ismip6_2dThermalForcing at monthly intervals), adds an ismip7_params stream reading ismip6shelfMelt_deltaT/_basin/_gamma0 from melt_params_path (mirroring the AIS convention, and consuming the output of the new build_3d_thermal_forcing step in landice/ismip7_forcing/ocean_thermal), and sets config_use_3d_thermal_forcing_for_face_melt = .true. melt_params_path is now validated as required when use_3d_thermal_forcing is true. streams.landice.template gates the TF stream's variable list and the new ismip7_params stream with Jinja2 conditionals on use_3d_thermal_forcing; verified both branches render as valid XML with the expected variables. Updates user and developer docs.
…SON config - _find_variable() returned dataset.variables[name] (a low-level xr.Variable), not a full DataArray. Variable.isel() doesn't accept drop=True, so passing a positional indexer dict together with drop=True raised 'cannot specify both keyword and positional arguments to .isel'. Return dataset[name] instead. - calibrate_regional_delta_t raised when a region had zero initially-floating cells, aborting the whole calibration. Skip such regions instead: print a diagnostic, leave their monthly means as NaN, set deltaT=0, and leave achieved melt as NaN (undefined, since there's no floating ice to melt). Also add greenland_3d_tf_config.json, a real example JSON config for the build_3d_thermal_forcing step (region-mask/EN4/GeoJSON paths only; mesh, forcing_2d, output, and diagnostics are always injected by the step), and wire it into ismip7_forcing_ocx_gis.cfg via process_ocean_thermal_3d and [ismip7_ocean_thermal_3d] config_file.
scipy's classic/CDF-2 (NETCDF3_64BIT) writer produces a header that netCDF-C (ncdump, MALI) rejects with "Unknown file format" when a dataset mixes a record (unlimited-Time) variable with a scalar (0-D) variable, as in this output (ismip6shelfMelt_gamma0 alongside xtime and ismip6shelfMelt_3dThermalForcing). scipy can read its own broken file back, masking the problem. Switch to engine="netcdf4", which writes a conformant CDF-2 file and still streams record-by-record.
deltaT/gamma0 are calibrated once against OCX and must be held fixed for every ESM scenario (recalibrating per-ESM would force all their historical mean melt rates toward the same targets and erase the differences between ESM ocean forcings). Config gains scenario and melt_params_file; Config.calibrate_delta_t gates calibration to scenario == OCX. Non-OCX runs load deltaT/gamma0/basin from an existing melt_params_file instead of recalibrating, validating that the stored basin layout and gamma0 match. Also split ismip6shelfMelt_basin/gamma0/deltaT out of the time-varying 3-D thermal forcing output into their own file (write_melt_params/ read_melt_params), independent of ismip6shelfMelt_3dThermalForcing. This matches how ismip7_run already expects a separate melt_params_path file (see ismip7_gris streams.landice.template ismip7_params stream). build_3d_thermal_forcing.py always points melt_params_file at the OCX output directory, regardless of which scenario is currently being processed.
The xtime char array from the 2-D forcing was carried over and run through xarray's character coder a second time, splitting each byte into a spurious length-1 dimension. Build xtime fresh from the decoded timestamps following the existing compass idiom (ljust(64) strings -> dtype 'S') and let xarray encode the StrLen char dimension, so xtime is written as a proper char(Time, StrLen) variable.
A 285-year monthly 3D forcing field would be ~113 GB in a single file. Write the output in blocks of output_years_per_file years (default 10), one file per block named by its block start year (..._<YYYY>.nc), which the run side addresses with a $Y filename_template plus a filename_interval. Also correct the GrIS 3D thermal-forcing cadence to monthly on the run side; the annual filename cadence only applies to Antarctica, whose ISMIP7 3D forcing is provided annually. - greenland_3d.py: add Config.output_years_per_file; year_chunks and chunk_output_path helpers; write one file per block via _write_forcing_chunk. - greenland_3d_tf_config.json: add output_years_per_file example. - ismip7_gris/set_up_experiment.py: symlink the whole chunk series, derive filename_interval and reference_time, force monthly TF. - ismip7_gris/streams.landice.template: templated TF filename_interval and reference_time.
Re-running the OCX step failed hard because the melt-params file (and forcing chunks) already existed with overwrite=false. Treat overwrite= false as "reuse existing, skip" so compass re-runs are idempotent and resumable (each chunk is written atomically via a .partial rename, so a present file is complete). overwrite=true still forces regeneration.
Writing NETCDF3_64BIT directly with an unlimited Time dimension and multiple record variables is several times slower because the classic writer interleaves each record and writes with a stride (benchmarked at ~143 s vs ~26 s for a 4 GB, 10-year monthly chunk). Write HDF5 first (variables stored contiguously) then convert to the CDF-2 format MALI/pnetcdf requires with nccopy, which preserves the unlimited Time dimension. The HDF5 scratch file is removed after conversion.
Where the seafloor is deeper than source_max_depth_m the anchor is clamped, so TF_3d matches TF_2d at source_max_depth_m rather than the true seafloor. If the regional profile still slopes at that depth, the per-cell offset (and thus the whole reconstructed column) is biased. Emit a runtime warning reporting the count, marine-ice area fraction, and deepest seafloor, and distinguish regions that still slope at the bottom (biased) from those that are flat there (benign), pointing to raising ocean_vertical_grid.bottom_m and source_max_depth_m together.
Strong upper-ocean seasonality makes a single annual vertical profile misleading, so group the EN4 monthly profiles into seasons_per_year equal calendar blocks (default 4: Jan-Mar, Apr-Jun, Jul-Sep, Oct-Dec) and build one profile per season per region. Each output month uses the column shape of its season, anchored to that month's 2-D forcing at the seafloor (TF_3d = TF_2d at the anchor is preserved). Calibration and the diagnostics also use the per-month seasonal profile; the profile plots now draw the seasonal means (and season-colored monthly curves) instead of one annual mean, and regional_profiles.csv gains a season column.
mapped_lons is wrapped to [-180, 180] but mesh.lon_deg is in [0, 360], so the EN4 points and the mesh landed on opposite sides of the region-assignment map (the assignment itself is correct; the KD-tree uses convention-independent xyz). Wrap the mesh longitudes to [-180, 180] in the region-assignment and representative-field maps so the two overlay.
Changes gamma0/deltaT calibration from single-stage (calibrate deltaT with fixed gamma0) to two-stage: Stage 1: Calibrate global gamma0 with all deltaT=0 - Uses only regions with non-None target melt rates - Minimizes mean squared error across active regions Stage 2 (optional): Calibrate deltaT per region with fixed gamma0 - Controlled by new calibrate_deltaT flag (defaults to false) - Only calibrates regions with non-None targets - Regions without targets keep deltaT=0 Key changes: - Add calibrate_deltaT boolean to Config (in calibration section of JSON) - Support None/null values in regional_melt_targets_m_per_yr - Add calibrate_gamma0() function using scipy.optimize.minimize_scalar - Refactor calibration into compute_regional_monthly_means() and calibrate_parameters() - Update write_melt_params() to accept calibrated gamma0 parameter - Update read_melt_params() to return (gamma0, deltaT) tuple - Enhanced reporting shows calibrated gamma0 and achieved melt for all regions regardless of calibration participation - Update example config to demonstrate None targets and calibrate_deltaT flag Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
The Python 3.13 CI build was failing because Sphinx couldn't fetch the xarray intersphinx inventory from the obsolete http://xarray.pydata.org URL, causing a connection reset error. With the -W flag (warnings as errors) in html-strict target, this terminated the build. Updated all intersphinx mappings to use HTTPS and the current xarray documentation domain (docs.xarray.dev). Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Replace manual get_params/forcing_group/source logic with resolve_ocean_source, so the step uses the same label and forcing_group that ProcessThermalForcing wrote. The 2D/3D filename contract now matches what process_thermal_forcing produces (mesh_<2d/3d>ThermalForcing_<label>_). This is a semantic-conflict resolution: build_3d originally mirrored pre-MPAS-Dev#997 source logic; now it consumes the MPAS-Dev#997 ForcingSource object. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
- Explicitly check for _YYYY.nc pattern match before calling .group(1), exiting with a clear error instead of raising AttributeError on None. - Validate that all adjacent chunk pairs have the same year spacing, not just the first two, before deriving tf_filename_interval. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
GrIS 3D thermal forcing is read at monthly intervals (matching 2D), with multiple year-block chunk files addressed via a $Y filename template. The annual cadence applies only to Antarctica. This was corrected in MPAS-Dev#979 commit 5b85afa but the doc wasn't fully updated during the rebase. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Both 2D and 3D GrIS thermal forcing use monthly input_interval. The 3D forcing is chunked across year-block files addressed with a $Y template and filename_interval. Only AIS uses annual thermal forcing intervals. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
trhille
force-pushed
the
landice/create_3d_ismip7_greenland_tf
branch
from
September 25, 2026 22:13
ff9e3c3 to
6914ad3
Compare
Contributor
There was a problem hiding this comment.
Copilot review overview
🟡 Changes recommended
Calibration behavior, basin handling, dependencies, diagnostics, documentation, and test coverage have unresolved issues.
Get a fresh assessment by requesting another Copilot review.
Review effort: Balanced
Findings: 2
Open (11)
What changed in this PR
Adds optional 3-D Greenland thermal-forcing generation from ISMIP7 2-D forcing and EN4 profiles, plus MALI runtime integration.
Changes:
- Adds regional profile construction, calibration, diagnostics, and chunked NetCDF output.
- Adds optional GrIS 3-D forcing streams and runtime configuration.
- Updates forcing filenames and documentation.
| File | Description |
|---|---|
docs/users_guide/landice/test_groups/ismip7_run.rst |
Documents GrIS 3-D runtime forcing. |
docs/users_guide/landice/test_groups/ismip7_forcing.rst |
Documents 3-D forcing generation. |
docs/developers_guide/landice/test_groups/ismip7_run.rst |
Describes runtime implementation. |
docs/developers_guide/landice/test_groups/ismip7_forcing.rst |
Describes processing internals. |
docs/developers_guide/landice/api.rst |
Adds new API entries. |
docs/conf.py |
Updates intersphinx URLs. |
compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template |
Adds conditional 3-D streams. |
compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py |
Configures chunked 3-D forcing. |
compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg |
Adds the 3-D forcing option. |
compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg |
Adds test configuration defaults. |
compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py |
Updates AIS filename matching. |
compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py |
Distinguishes 2-D and 3-D filenames. |
compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py |
Implements the new processing pipeline. |
compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py |
Integrates processing as a Compass step. |
compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py |
Registers the new step. |
compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg |
Adds shared 3-D options. |
compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg |
Adds OCX Greenland configuration. |
compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json |
Supplies EN4 and calibration parameters. |
💡 Configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| f"{unassigned.size} cells have no region; first indices: " | ||
| f"{unassigned[:10].tolist()}" | ||
| ) | ||
| return np.argmax(membership, axis=1).astype(np.int32) + 1 |
Comment on lines
+1369
to
+1373
| # Stage 1: Calibrate gamma0 with deltaT=0 using only non-NaN targets | ||
| print("\n=== Stage 1: Calibrating global gamma0 (with deltaT=0) ===") | ||
| gamma0 = calibrate_gamma0( | ||
| monthly_means, cfg.regional_melt_targets_m_per_yr, coefficient | ||
| ) |
Comment on lines
+2316
to
+2319
| def run(cfg: Config, logger, prepare_only: bool = False) -> None: | ||
| mesh, basin_ids = load_mesh_and_basins(cfg) | ||
| profiles = build_regional_profiles(cfg, mesh, basin_ids, logger) | ||
| warn_if_seafloor_below_max_depth(cfg, mesh, basin_ids, profiles, logger) |
Comment on lines
+525
to
+526
| ocean_thermal.greenland_3d.Config | ||
| ocean_thermal.greenland_3d.run |
Comment on lines
+157
to
+158
| ``ismip6shelfMelt_basin``. The multi-gigabyte field is streamed record by | ||
| record (dask, ``scipy`` engine, ``NETCDF3_64BIT``, ``.partial``-then-rename). |
Comment on lines
+107
to
+108
| {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_2dThermalForcing_{source}_{scenario}_{years}.nc (GrIS) | ||
| {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_3dThermalForcing_{source}_{scenario}_{years}.nc (AIS; optional GrIS 3-D) |
Comment on lines
+372
to
+373
| ``process_ocean_thermal_3d = true``) converts the GrIS 2D thermal forcing into | ||
| a 30-level 3D field for MALI's nonlocal (Jourdain et al. 2020) melt scheme. It |
Comment on lines
+250
to
+252
| When ``use_3d_thermal_forcing`` is ``false`` (the default), GrIS looks for | ||
| ``*2dThermalForcing_*.nc`` and reads ``ismip6_2dThermalForcing`` at monthly | ||
| intervals. When ``true``, it instead looks for ``*3dThermalForcing_*.nc`` and |
Comment on lines
+277
to
+278
| * ``ismip6shelfMelt_3dThermalForcing`` (AIS; GrIS when | ||
| ``use_3d_thermal_forcing = true``) |
Add gis_contShelfExtent_EPSG4326.geojson defining the Greenland continental shelf extent in WGS84 coordinates. This is now the default en4_source_region_geojson, eliminating the need to specify it in the JSON config unless users want to override with a custom region. The GeoJSON contains a single Polygon covering the Greenland continental shelf, derived from GSFC drainage basin outlines, used to filter EN4 source points when building regional thermal forcing profiles. Updated greenland_3d_tf_config.json to remove the en4_source_region_geojson path (now optional with bundled default) and added a comment explaining the default behavior. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
The gamma0_m_per_yr field in the calibration config appeared to set the Jourdain nonlocal-scheme gamma0 but was not actually used during calibration. calibrate_gamma0() computes gamma0 via bounded optimization (100-100000 m/yr) and never reads cfg.gamma0_m_per_yr. The field was only used for: 1. A positivity validation check (now removed) 2. Being written to the deltaT_calibration.json diagnostic (misleading - it reported the config value 14500.0 while the melt-params file got the actually-calibrated value) ESM (non-OCX) runs read gamma0 from the melt-params file via read_melt_params and never used this config field either. Changes: - Remove gamma0_m_per_yr field from Config dataclass - Remove parsing from from_json() - Remove positivity validation check - Fix diagnostic: write_diagnostics() now reports the calibrated gamma0 (which it already receives as a parameter) instead of cfg.gamma0_m_per_yr - Remove gamma0_m_per_yr from example JSON config After this change, deltaT_calibration.json will correctly report the gamma0 that was actually calibrated and written to the melt-params file. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Changed GrIS run setup so ismip6_2dThermalForcing is ALWAYS read (the authoritative ISMIP7 2D forcing), with ismip6shelfMelt_3dThermalForcing added alongside when use_3d_thermal_forcing = true, instead of the previous either/or behavior. Rationale: The 2D forcing is the true ISMIP7 forcing; the 3D product is a separate derived product and should not replace it. MALI consumes both fields simultaneously. Changes: streams.landice.template: - ismip7_TF stream now unconditionally reads ismip6_2dThermalForcing from a single file with static reference_time=2000-01-01_00:00:00 - New ismip7_TF_3d stream (gated on use_3d_thermal_forcing) reads ismip6shelfMelt_3dThermalForcing + ismip6shelfMelt_zOcean from chunked $Y series with filename_interval set_up_experiment.py: - Always glob *2dThermalForcing_*.nc and symlink for 2D stream - When use_3d_thermal_forcing=true, additionally glob *3dThermalForcing_*.nc, symlink the chunk series, derive tf_3d_filename_interval from uniform chunk spacing - CTRL: add ctrl_tf_3d_climatology_path support, error if use_3d but path NotAvailable - stream_replacements: add input_file_TF_3d_forcing and tf_3d_filename_interval; drop old tf_reference_time/tf_filename_interval (now hardcoded static in template) Config files (ismip7_gris.cfg, ismip7_gris_test.cfg): - Add ctrl_tf_3d_climatology_path option (default NotAvailable), needed for CTRL + use_3d_thermal_forcing Docs (users_guide and developers_guide ismip7_run.rst): - Correct to say 2D is always read; 3D is read additionally when flag is true - Update Forcing Streams section to show 2D always present, 3D monthly when enabled - Clarify both forcings delivered to MALI simultaneously Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Addressed review comments MPAS-Dev#1, MPAS-Dev#2, MPAS-Dev#3, MPAS-Dev#4, MPAS-Dev#5, MPAS-Dev#6, MPAS-Dev#7, MPAS-Dev#8 from the Copilot review. Deferred MPAS-Dev#9 (Dask dependency) and MPAS-Dev#10 (test suite) as low priority. CODE FIXES: region) to the nearest region by great-circle distance (cKDTree), instead of silently assigning them to basin 1 via argmax. Overlapping cells (multiple regions) still pick the first matching region. Added lat/lon parameters to build_basin_ids and updated the call site in load_mesh_and_basins. greenland_3d.py so only the public API (Config and run) remains non- underscored. This satisfies the repository convention that non-underscore module-scope symbols should be documented in api.rst. ice) with a finite JSON target would crash calibrate_gamma0 instead of being excluded. Now NaN out empty regions in the target array before calibration, so active_regions correctly excludes both NaN targets and empty regions. DOC FIXES: (it uses forcing_interval_annual). GrIS 3D stays in Monthly (correct). {mesh}_2dThermalForcing_*.nc (GrIS) and {mesh}_3dThermalForcing_*.nc (AIS; optional GrIS), and fixed lowercase smb→SMB. locations (user/dev ismip7_forcing.rst) since the checked-in JSON requests 20 levels and it's configurable. NETCDF3_64BIT" to "netCDF4 engine → nccopy -k nc6 for CDF-2" (actual implementation uses netCDF4 + nccopy conversion). block (required by ismip7_run as melt_params_path but was undocumented). Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Collaborator
Author
TestingI ran this with the following .cfg file: click to see .cfgand this json file: click to see .jsonHere is the log output written to click to see log fileAnd here are the outputs created in the |
Remove code that overwrote the default namelist setting of config_use_3d_thermal_forcing_for_face_melt = .false. We always want to use the 2d thermal forcing provided by ISMIP7 for face-melting, since constructing the 3d forcing has the potential to introduce errors.
Use mean absolute error instead of mean squared error to calibrate gamma0 to desired melt rates. Using MSE led to very low melt rates in the sector with the second-largest amount of floating ice.
Replace optimizer with closed-form solve: gamma0 calibrated so the area-weighted mean melt across finite-target regions equals the area-weighted mean of their targets. Fixes two bugs in the old L1 fit: no across-region area weighting (85 km² sliver = 2551 km² shelf) and median-collapse onto one region (north_east pinned at 10.00000 exactly). Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Collaborator
Author
|
After 997667e, I re-ran the calibration using the same .cfg settings but a modified json file: The calibration results were: Diagnostic plots: |
Member
|
I plotted up profiles of TF in the generated .nc files. I plot 5 random ocean cells in each of the GIS regions: Profiles for 2010-01 time level: 
|
Member
|
And here are averaged profiles for each region for the 4 seasons (indices 0, 3, 6, 9 of 
|
matthewhoffman
approved these changes
Sep 27, 2026
matthewhoffman
left a comment
matthewhoffman
left a comment
Member
There was a problem hiding this comment.
Approving based on demonstrating that the profiles are doing what they are supposed to!
Fix a small discrepancy between the freezing-point slope used in the gamma0 calculation and the value used in MALI. Also add one deep ocean level to match bottom of observational data set.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Create 3d Greenland thermal forcing from ISMIP7-provided 2d and EN4 data. Note that this is currently branched from #978. We will rebase once that is merged.
Checklist
api.rst) has any new or modified class, method and/or functions listedTestingin this PR) any testing that was used to verify the changes