diff --git a/CHANGELOG.md b/CHANGELOG.md index 0bdeaddc..de3ef997 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ SPDX-License-Identifier: CC-BY-4.0 # develop +* Fixed off-by-one error in grid shifts that can lead to inaccuracies for long-running simulations. https://github.com/EBFMorg/EBFM/pull/164 * Add a 10-km test case for Greenland forced with CARRA2 meteorological data for May 1990 (`examples/Greenland_CARRA2`). https://github.com/EBFMorg/EBFM/pull/128 * Add feature to use NetCDF meshes with DEM as input (`--netcdf-mesh`) * Add option to set missing elevations of NetCDF meshes to 0 m (`--fill-missing-elevation`) diff --git a/src/ebfm/core/LOOP_SNOW.py b/src/ebfm/core/LOOP_SNOW.py index 21c07416..1b34e80f 100644 --- a/src/ebfm/core/LOOP_SNOW.py +++ b/src/ebfm/core/LOOP_SNOW.py @@ -147,10 +147,12 @@ def snowfall_and_deposition(): subD_old = OUT["subD"][idx_shift] subW_old = OUT["subW"][idx_shift] - OUT["subZ"][idx_shift, 2 : nl - 1] = subZ_old[:, 1 : nl - 2] - OUT["subT"][idx_shift, 2 : nl - 1] = subT_old[:, 1 : nl - 2] - OUT["subD"][idx_shift, 2 : nl - 1] = subD_old[:, 1 : nl - 2] - OUT["subW"][idx_shift, 2 : nl - 1] = subW_old[:, 1 : nl - 2] + # Shift the column down by one: the layers at indices 1..nl-2 move to + # indices 2..nl-1. The layer at index nl-1 falls out of the domain. + OUT["subZ"][idx_shift, 2:nl] = subZ_old[:, 1 : nl - 1] + OUT["subT"][idx_shift, 2:nl] = subT_old[:, 1 : nl - 1] + OUT["subD"][idx_shift, 2:nl] = subD_old[:, 1 : nl - 1] + OUT["subW"][idx_shift, 2:nl] = subW_old[:, 1 : nl - 1] OUT["subZ"][idx_shift, 1] = max_subZ OUT["subZ"][idx_shift, 0] = (subZ_old[:, 0] + shift_amount[idx_shift]) - max_subZ @@ -168,7 +170,8 @@ def snowfall_and_deposition(): OUT["subW"][idx_shift, 1] = subW_old[:, 0] OUT["subW"][idx_shift, 0] = 0.0 - # Update runoff for shifted layers + # The irreducible water in the previously deepest layer (now removed) is added to runoff_irr_deep for + # mass conservation. OUT["runoff_irr_deep"][idx_shift] += subW_old[:, nl - 1] return _SUCCESS @@ -238,10 +241,12 @@ def melt_sublimation(): subD_old = OUT["subD"][idx_shift] subW_old = OUT["subW"][idx_shift] - OUT["subZ"][idx_shift, 1 : nl - 2] = subZ_old[:, 2 : nl - 1] - OUT["subT"][idx_shift, 1 : nl - 2] = subT_old[:, 2 : nl - 1] - OUT["subD"][idx_shift, 1 : nl - 2] = subD_old[:, 2 : nl - 1] - OUT["subW"][idx_shift, 1 : nl - 2] = subW_old[:, 2 : nl - 1] + # Shift the column up by one: the layers at indices 2..nl-1 move to + # indices 1..nl-2. + OUT["subZ"][idx_shift, 1 : nl - 1] = subZ_old[:, 2:nl] + OUT["subT"][idx_shift, 1 : nl - 1] = subT_old[:, 2:nl] + OUT["subD"][idx_shift, 1 : nl - 1] = subD_old[:, 2:nl] + OUT["subW"][idx_shift, 1 : nl - 1] = subW_old[:, 2:nl] subZ_top_new = subZ_old[:, 0] + subZ_old[:, 1] + shift_amount[idx_shift] OUT["subZ"][idx_shift, 0] = subZ_top_new @@ -249,10 +254,11 @@ def melt_sublimation(): OUT["subD"][idx_shift, 0] = subD_old[:, 1] # Liquid water scales with the remaining fraction of the former layer 1 OUT["subW"][idx_shift, 0] = subW_old[:, 1] * (subZ_top_new / subZ_old[:, 1]) + + # Index nl-1 becomes a fresh layer at the base: it inherits the temperature and density of the + # pre-shift deepest layer, starts dry, and gets the full deep-layer thickness. OUT["subT"][idx_shift, nl - 1] = subT_old[:, nl - 1] OUT["subW"][idx_shift, nl - 1] = 0.0 - - # Update the deepest layer properties (vectorized over shift rows) if column.doubledepth: OUT["subZ"][idx_shift, nl - 1] = 2.0 ** len(column.split) * column.max_subZ else: diff --git a/tests/core/test_loop_snow_layer_shift.py b/tests/core/test_loop_snow_layer_shift.py new file mode 100644 index 00000000..88066053 --- /dev/null +++ b/tests/core/test_loop_snow_layer_shift.py @@ -0,0 +1,256 @@ +# SPDX-FileCopyrightText: 2026 EBFM Authors +# +# SPDX-License-Identifier: BSD-3-Clause +# +# Content partially generated with the assistance of AI tools. +# Claude Code: Opus 5 + +"""The vertical grid shifts in LOOP_SNOW. + +Two blocks move the subsurface column by one layer index: snowfall_and_deposition +pushes it down when a new top layer is started, and melt_sublimation pulls it up +when the top layer is gone. Both must move *every* interior layer, and the layer +that drops off the bottom is the one snowfall_and_deposition credits to +runoff_irr_deep as water leaving the domain. + +These tests use a column whose per-layer values are all distinct, and in +particular whose two deepest layers differ. That is essential: the model's own +initial state gives every deep layer the same thickness, temperature, density and +water content (INIT sets subD = C["Dice"] throughout and subZ[34:] = 0.8), and +duplicating or dropping a layer identical to its neighbour is numerically +invisible. A test built on a uniform deep column passes whether the shifts are +right or wrong. + +subW is the tracer of choice, because in this setup nothing but the shift can +touch it: the column sits exactly at the melting point so it has no cold content +and no water can refreeze, no rain or melt adds water, and every layer's water +stays below the mliqmax cap that compaction applies. That makes the assertions on +subW exact rather than approximate. +""" + +import unittest + +import numpy as np + +from ebfm.core.config import ColumnDiscretizationConfig + +GPSUM, NL = 2, 50 +# Layer thickness (m). The top layer is a full max_subZ, so any snowfall starts a +# new one and forces a downward shift. +MAX_SUBZ = 0.1 +THICKNESS = 0.5 +# Distinct per layer, and inside [C["Dfirn"], C["Dice"]) so only firn compaction +# applies. The 3.0 kg m-3 spacing is what makes a misplaced layer detectable. +DENSITY_STEP = 3.0 +DENSITY_BASE = 600.0 +# Distinct per layer, and far below the mliqmax of these layers (9.6 kg m-2 at the +# base rising to 13.5 at the top), so compaction's cap leaves them untouched. +WATER_BASE, WATER_STEP = 1.0, 0.1 + + +def _make_case(C): + """A column set at the melting point whose every layer differs from its neighbours.""" + subZ = np.full(NL, THICKNESS) + subZ[0] = MAX_SUBZ + subD = DENSITY_BASE + DENSITY_STEP * np.arange(NL) + subW = WATER_BASE + WATER_STEP * np.arange(NL) + + grid = { + "gpsum": GPSUM, + "nl": NL, + "max_subZ": MAX_SUBZ, + "doubledepth": False, + "split": np.array([15]), + # Layer merging/splitting does its own shifts and its own runoff_irr_deep + # credit; switch it off so these tests see the snowfall/melt shifts alone. + "mask": np.zeros(GPSUM, dtype=int), + } + OUT = { + "subZ": np.tile(subZ, (GPSUM, 1)), + "subD": np.tile(subD, (GPSUM, 1)), + "subW": np.tile(subW, (GPSUM, 1)), + "subT": np.full((GPSUM, NL), C["T0"]), + "subTmean": np.full((GPSUM, NL), C["T0"]), + "subS": np.zeros((GPSUM, NL)), + "surfH": np.zeros(GPSUM), + "Tsurf": np.full(GPSUM, C["T0"]), + "melt": np.zeros(GPSUM), + "moist_deposition": np.zeros(GPSUM), + "moist_sublimation": np.zeros(GPSUM), + "moist_condensation": np.zeros(GPSUM), + "moist_evaporation": np.zeros(GPSUM), + "runoff_irr_deep_mean": np.zeros(GPSUM), + } + yearsnow = np.full(GPSUM, 500.0) + IN = { + "T": np.full(GPSUM, C["T0"]), + "WS": np.full(GPSUM, 4.0), + "snow": np.zeros(GPSUM), + "rain": np.zeros(GPSUM), + "yearsnow": np.tile(yearsnow[:, None], (1, NL)), + "logyearsnow": np.tile(np.log(yearsnow)[:, None], (1, NL)), + } + return grid, OUT, IN + + +class TestColumnShift(unittest.TestCase): + def setUp(self): + from ebfm.core import INIT, LOOP_SNOW + + self.LOOP_SNOW = LOOP_SNOW + self.C = INIT.init_constants() + self.phys = {"snow_compaction": "firn+snow", "percolation": "normal"} + self.column = ColumnDiscretizationConfig() + + def _step(self, grid, OUT, IN): + self.LOOP_SNOW.main(self.C, OUT, IN, 1.0 / 24.0, grid, self.phys, self.column) + + def _snowfall_case(self): + """A case that performs exactly one downward shift. + + The fresh snow is sized from the fresh-snow density the model will compute, so + that it overfills the top layer without filling a second one: one pass of the + shift loop, whose result is a plain permutation of the column. + """ + grid, OUT, IN = _make_case(self.C) + Dfresh = 50 + 1.7 * 17 ** (3 / 2) + 266.86 * (0.5 * (1 + np.tanh(4.0 / 5))) ** 8.8 + IN["snow"][:] = 0.9 * MAX_SUBZ * Dfresh / self.C["Dwater"] + return grid, OUT, IN + + def _sublimation_case(self): + """A case that performs exactly one upward shift. + + Sublimation, not melt, removes the top layer: it takes away exactly as much mass + as the top layer holds, and unlike melt it produces no liquid water, so subW + stays an exact tracer. Both drive the same block through the same `shift_tot`. + """ + grid, OUT, IN = _make_case(self.C) + OUT["moist_sublimation"][:] = OUT["subD"][:, 0] * OUT["subZ"][:, 0] * 1e-3 + return grid, OUT, IN + + def test_snowfall_shifts_every_layer_down(self): + """Every layer moves down one index; the deepest layer leaves the domain. + + The layer that must not survive is the deepest one: if it is still there + afterwards, the column kept a layer that the runoff bookkeeping already + counted as gone. + """ + grid, OUT, IN = self._snowfall_case() + subW_old = OUT["subW"].copy() + + self._step(grid, OUT, IN) + + want = np.empty((GPSUM, NL)) + want[:, 0] = 0.0 # new top layer is fresh snow, no water + want[:, 1] = subW_old[:, 0] # former top layer, handed down whole + want[:, 2:NL] = subW_old[:, 1 : NL - 1] # everything else moves down one index + np.testing.assert_array_equal( + OUT["subW"], + want, + err_msg="subW is not the column shifted down by one index", + ) + + def test_sublimation_shifts_every_layer_up(self): + """Every layer moves up one index; a fresh, dry layer appears at the base.""" + grid, OUT, IN = self._sublimation_case() + subW_old = OUT["subW"].copy() + + self._step(grid, OUT, IN) + + want = np.empty((GPSUM, NL)) + want[:, 0] = subW_old[:, 1] # former layer 1 becomes the whole top layer + want[:, 1 : NL - 1] = subW_old[:, 2:NL] # everything else moves up one index + want[:, NL - 1] = 0.0 # fresh layer added at the base + np.testing.assert_allclose( + OUT["subW"], + want, + rtol=1e-12, + atol=0.0, + err_msg="subW is not the column shifted up by one index", + ) + + def test_snowfall_shifts_density_down(self): + """The density profile rides the same shift as the water profile. + + Firn compaction perturbs subD by ~1e-3 kg m-3 over one hourly step, three + orders of magnitude below the layer-to-layer spacing, so a misplaced layer is + unambiguous at this tolerance. + """ + grid, OUT, IN = self._snowfall_case() + subD_old = OUT["subD"].copy() + + self._step(grid, OUT, IN) + + np.testing.assert_allclose( + OUT["subD"][:, 2:NL], + subD_old[:, 1 : NL - 1], + rtol=1e-4, + err_msg="subD is not the column shifted down by one index", + ) + + def test_sublimation_shifts_density_up(self): + """The density profile rides the same shift as the water profile.""" + grid, OUT, IN = self._sublimation_case() + subD_old = OUT["subD"].copy() + + self._step(grid, OUT, IN) + + np.testing.assert_allclose( + OUT["subD"][:, 1 : NL - 1], + subD_old[:, 2:NL], + rtol=1e-4, + err_msg="subD is not the column shifted up by one index", + ) + + def test_snowfall_shift_conserves_column_mass(self): + """Column mass changes by exactly the snow added minus the layer that left. + + Nothing else in this configuration moves mass: compaction preserves subD*subZ + by construction, heat conduction touches neither, and with no rain, no melt and + no cold content there is no refreezing. So the whole residual is the shift's. + """ + grid, OUT, IN = self._snowfall_case() + mass_before = np.sum(OUT["subD"] * OUT["subZ"], axis=1) + # The fresh snow mass is the snowfall in water equivalent, whatever density the + # model gives it: shift_snowfall * Dfreshsnow == snow * Dwater. + mass_added = IN["snow"] * self.C["Dwater"] + mass_left = OUT["subD"][:, NL - 1] * OUT["subZ"][:, NL - 1] + + self._step(grid, OUT, IN) + + mass_after = np.sum(OUT["subD"] * OUT["subZ"], axis=1) + np.testing.assert_allclose( + mass_after - mass_before, + mass_added - mass_left, + atol=1e-6, + err_msg="column mass across the shift does not balance snowfall in minus deepest layer out", + ) + + def test_water_credited_to_deep_runoff_actually_leaves(self): + """runoff_irr_deep must be credited the water that really left the column. + + This is the bookkeeping the shift has to agree with: the block credits the + deepest layer's water as having left the domain, so the deepest layer's water + is what the shift has to discard. + """ + grid, OUT, IN = self._snowfall_case() + water_before = np.sum(OUT["subW"], axis=1) + + self._step(grid, OUT, IN) + + water_after = np.sum(OUT["subW"], axis=1) + # runoff() has already rescaled runoff_irr_deep into the running mean, so read + # the credit back off the mean, which started at zero. + weight = (1.0 / 24.0) / self.C["yeardays"] + credited = OUT["runoff_irr_deep_mean"] / weight + + np.testing.assert_allclose( + credited, + water_before - water_after, + rtol=1e-12, + err_msg="water credited to runoff_irr_deep differs from the water that left the column", + ) + + +if __name__ == "__main__": + unittest.main()