From fe1635be653d1d2b5194b0925943c7a0318e2c9c Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Mon, 10 Nov 2025 18:33:29 -0500
Subject: [PATCH 1/3] Store OGGM flowline diagnostics

---
 .../postproc_binned_subannual_thick.py        | 48 ++++++-----
 pygem/bin/run/run_calibration.py              | 38 ++++++---
 pygem/bin/run/run_simulation.py               | 79 +++++++++++--------
 pygem/massbalance.py                          |  7 ++
 pygem/output.py                               | 26 ++++++
 5 files changed, 125 insertions(+), 73 deletions(-)

diff --git a/pygem/bin/postproc/postproc_binned_subannual_thick.py b/pygem/bin/postproc/postproc_binned_subannual_thick.py
index 13479ae3..e917776a 100644
--- a/pygem/bin/postproc/postproc_binned_subannual_thick.py
+++ b/pygem/bin/postproc/postproc_binned_subannual_thick.py
@@ -67,7 +67,13 @@ def getparser():
 
 
 def get_binned_subannual(
-    bin_massbalclim, bin_mass_annual, bin_thick_annual, dates_subannual, dates_annual, debug=False
+    bin_massbalclim,
+    bin_mass_annual,
+    bin_thick_annual,
+    bin_flux_divergence_annual,
+    dates_subannual,
+    dates_annual,
+    debug=False,
 ):
     """
     funciton to calculate the subannual binned ice thickness and mass
@@ -79,10 +85,8 @@ def get_binned_subannual(
     thus, subannual thickness and mass is determined assuming
     the flux divergence is constant throughout the year.
 
-    annual flux divergence is first estimated by combining the annual binned change in ice
-    thickness and the annual binned mass balance. then, assume flux divergence is constant
-    throughout the year (divide annual by the number of steps in the binned climatic mass
-    balance to get subannual flux divergence).
+    assumes flux divergence is constant throughout the year, dividing annual by the number of steps
+    in the binned climatic mass balance to get subannual flux divergence).
 
     subannual binned flux divergence can then be combined with
     subannual binned climatic mass balance to get subannual binned change in ice thickness and mass.
@@ -99,6 +103,9 @@ def get_binned_subannual(
     bin_thick_annual : ndarray
        annual binned glacier thickness [m ice]
         shape : [#glac, #elevbins, #years]
+    bin_flux_divergence_annual : ndarray
+       annual binned flux divergence [m ice]
+        shape : [#glac, #elevbins, #years]
     dates_subannual : array-like of datetime-like
         dates associated with `bin_massbalclim` (subannual)
     dates_annual : array-like of datetime-like
@@ -130,33 +137,23 @@ def get_binned_subannual(
     rho_i = pygem_prms['constants']['density_ice']
     bin_massbalclim_ice = bin_massbalclim * (rho_w / rho_i)
 
-    # --- Step 2: compute annual cumulative mass balance ---
-    # Initialize annual cumulative mass balance (exclude last year for flux calculation)
-    bin_massbalclim_annual = np.zeros((n_glac, n_bins, nyrs))
-    for i, year in enumerate(yrs):
-        idx = np.where(years_subannual == year)[0]
-        bin_massbalclim_annual[:, :, i] = bin_massbalclim_ice[:, :, idx].sum(axis=-1)
-
-    # --- Step 3: compute annual thickness change ---
-    bin_delta_thick_annual = np.diff(bin_thick_annual, axis=-1)  # [m ice yr^-1]
-
-    # --- Step 4: compute annual flux divergence ---
-    bin_flux_divergence_annual = bin_massbalclim_annual - bin_delta_thick_annual  # [m ice yr^-1]
-
-    # --- Step 5: expand flux divergence to subannual steps ---
+    # --- Step 2: expand flux divergence to subannual steps ---
+    # assume flux divergence is constant throughout the year
+    # (divide annual by the number of steps in the binned climatic mass balance to get subannual flux divergence)
+    bin_flux_divergence_annual = bin_flux_divergence_annual[:, :, 1:]
     bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice)
     for i, year in enumerate(yrs):
         idx = np.where(years_subannual == year)[0]
         bin_flux_divergence_subannual[:, :, idx] = bin_flux_divergence_annual[:, :, i][:, :, np.newaxis] / len(idx)
 
-    # --- Step 6: compute subannual thickness change ---
+    # --- Step 3: compute subannual thickness change ---
     bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual
 
-    # --- Step 7: calculate subannual thickness ---
-    running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1)
+    # --- Step 4: calculate subannual thickness ---
+    running_bin_delta_thick_subannual = np.nancumsum(bin_delta_thick_subannual, axis=-1)
     bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, :, 0][:, :, np.newaxis]
 
-    # --- Step 8: compute glacier volume and area on subannual timestep ---
+    # --- Step 5: compute glacier volume and area on subannual timestep ---
     bin_volume_annual = bin_mass_annual / rho_i  # annual volume [m^3] per bin
     bin_area_annual = np.divide(
         bin_volume_annual[:, :, 1:],  # exclude first year to match flux_div_annual
@@ -165,7 +162,7 @@ def get_binned_subannual(
         where=bin_thick_annual[:, :, 1:] > 0,
     )
 
-    # --- Step 9 : compute subannual glacier mass ---
+    # --- Step 6 : compute subannual glacier mass ---
     # First expand area to subannual steps
     bin_area_subannual = np.full(bin_massbalclim_ice.shape, np.nan)
     for i, year in enumerate(yrs):
@@ -175,7 +172,7 @@ def get_binned_subannual(
     # multiply by ice density to get subannual mass
     bin_mass_subannual = bin_thick_subannual * rho_i * bin_area_subannual
 
-    # --- Step 10: debug check ---
+    # --- Step 7: debug check ---
     if debug:
         for i, year in enumerate(yrs):
             # get last subannual index of that year
@@ -302,6 +299,7 @@ def run(simpath, debug):
             bin_massbalclim=binned_ds.bin_massbalclim.values,
             bin_mass_annual=binned_ds.bin_mass_annual.values,
             bin_thick_annual=binned_ds.bin_thick_annual.values,
+            bin_flux_divergence_annual=binned_ds.bin_flux_divergence_annual.values,
             dates_subannual=dates_subannual,
             dates_annual=dates_annual,
             debug=debug,
diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 1899e272..eef33cd4 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -317,6 +317,7 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
                 ev_model.mb_model.glac_wide_massbaltotal = (
                     ev_model.mb_model.glac_wide_massbaltotal + ev_model.mb_model.glac_wide_frontalablation
                 )
+        ev_model.mb_model.ensure_mass_conservation(diag)
     # safely catch any errors with dynamical run
     except Exception:
         ds = None
@@ -324,7 +325,7 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
     return mbmod, ds
 
 
-def calc_thick_change_1d(gdir, mbmod, ds):
+def calc_thick_change_1d(gdir, fls, mbmod, ds):
     """
     calculate binned change in ice thickness assuming constant annual flux divergence.
     sub-annual ice thickness is differenced at timesteps coincident with observations.
@@ -335,9 +336,6 @@ def calc_thick_change_1d(gdir, mbmod, ds):
     # grab components of interest
     bin_thick_annual = ds[0].thickness_m.values.T  # glacier thickness [m ice], (nbins, nyears)
 
-    # set any < 0 thickness to nan
-    bin_thick_annual[bin_thick_annual <= 0] = np.nan
-
     #  --- Step 1: convert mass balance from m w.e. to m ice  ---
     bin_massbalclim = mbmod.glac_bin_massbalclim  # climatic mass balance [m w.e.] per step
     # convert to m ice
@@ -358,12 +356,16 @@ def calc_thick_change_1d(gdir, mbmod, ds):
     bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual
 
     # --- Step 4: calculate subannual thickness = running thickness change + initial thickness---
-    running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1)
+    running_bin_delta_thick_subannual = np.nancumsum(bin_delta_thick_subannual, axis=-1)
     bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, 0][:, np.newaxis]
 
     # --- Step 5: rebin ---
     # get surface height at the specified reference year
-    ref_surface_height = ds[0].bed_h.values + ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values
+    ref_surface_height = (
+        fls[0].surface_h
+        - ds[0].thickness_m.isel(time=0).values
+        + ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values
+    )
     # aggregate model bin thicknesses as desired
     with warnings.catch_warnings():
         warnings.filterwarnings('ignore')
@@ -416,7 +418,7 @@ def mcmc_model_eval(
     if calib_elev_change_1d:
         mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls)
         # note, the binned thickness change is scaled by modeled density in mcmc.mbPosterior.log_likelihood() to calculate modeled surface elevation change
-        results['elev_change_1d'] = calc_thick_change_1d(gdir, mbmod, ds) if ds else float('-inf')
+        results['elev_change_1d'] = calc_thick_change_1d(gdir, fls, mbmod, ds) if ds else float('-inf')
 
     if mbfxn is not None:
         # grab current values from modelprms for the emulator
@@ -859,6 +861,19 @@ def run(list_packed_vars):
                         if not isinstance(gdir.elev_change_1d['dh_sigma'], int)
                         else gdir.elev_change_1d['dh_sigma']
                     )
+
+                    # optionally adjust ref_startyear based on earliest available elevation calibration data (must be <= 2000)
+                    if args.spinup:
+                        args.ref_startyear = min(
+                            2000, *(int(date[:4]) for pair in gdir.elev_change_1d['dates'] for date in pair)
+                        )
+                        # adjust dates table and climate data
+                        mask = gdir.dates_table['year'] >= args.ref_startyear
+                        gdir.dates_table = gdir.dates_table.loc[mask].reset_index(drop=True)
+                        nrows_rm = (~mask).sum()  # number of dropped rows
+                        for key in ('temp', 'tempstd', 'prec', 'lr'):
+                            gdir.historical_climate[key] = gdir.historical_climate[key][nrows_rm:]
+
                     # get observation period indices in model date_table
                     # create lookup dict (timestamp → index)
                     date_to_index = {d: i for i, d in enumerate(gdir.dates_table['date'])}
@@ -869,11 +884,6 @@ def run(list_packed_vars):
                         )
                         for start, end in gdir.elev_change_1d['dates']
                     ]
-                    # optionally adjust ref_startyear based on earliest available elevation calibration data (must be <= 2000)
-                    if args.spinup:
-                        args.ref_startyear = min(
-                            2000, *(int(date[:4]) for pair in gdir.elev_change_1d['dates'] for date in pair)
-                        )
 
                     # load calibrated calving_k values for tidewater glaciers
                     if gdir.is_tidewater and pygem_prms['setup']['include_frontalablation']:
@@ -1884,7 +1894,7 @@ def calc_mb_total_minelev(modelprms):
                     ]
                     # Melt
                     # energy available for melt [degC day]
-                    melt_energy_available = T_minelev * dates_table['days_in_step'].values
+                    melt_energy_available = T_minelev * gdir.dates_table['days_in_step'].values
                     melt_energy_available[melt_energy_available < 0] = 0
                     # assume all snow melt because anything more would melt underlying ice in lowermost bin
                     # SNOW MELT [m w.e.]
@@ -2135,6 +2145,8 @@ def rho_constraints(**kwargs):
                 )
                 # prepare export modelprms dictionary
                 modelprms_export = {}
+                # store reference period
+                modelprms_export['ref_period'] = (args.ref_startyear, args.ref_endyear)
                 # store model parameters and priors
                 modelprms_export['precgrad'] = [pygem_prms['sim']['params']['precgrad']]
                 modelprms_export['tsnow_threshold'] = [pygem_prms['sim']['params']['tsnow_threshold']]
diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py
index 3a3bc0fd..22ca9935 100755
--- a/pygem/bin/run/run_simulation.py
+++ b/pygem/bin/run/run_simulation.py
@@ -976,7 +976,7 @@ def run(list_packed_vars):
                                 nfls,
                                 y0=args.sim_startyear,
                                 mb_model=mbmod,
-                                glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
+                                glen_a=glen_a,
                                 fs=fs,
                                 is_tidewater=gdir.is_tidewater,
                                 water_level=water_level,
@@ -987,7 +987,7 @@ def run(list_packed_vars):
                                 nfls,
                                 y0=args.sim_startyear,
                                 mb_model=mbmod,
-                                glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
+                                glen_a=glen_a,
                                 fs=fs,
                             )
 
@@ -997,7 +997,7 @@ def run(list_packed_vars):
                                 ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}'
                             )
 
-                        diag = ev_model.run_until_and_store(args.sim_endyear + 1)
+                        diag, ds = ev_model.run_until_and_store(args.sim_endyear + 1, fl_diag_path=True)
                         ev_model.mb_model.glac_wide_volume_annual[-1] = diag.volume_m3[-1]
                         ev_model.mb_model.glac_wide_area_annual[-1] = diag.area_m2[-1]
 
@@ -1231,8 +1231,19 @@ def run(list_packed_vars):
                         output_offglac_melt_steps[:, n_iter] = mbmod.offglac_wide_melt
                         output_offglac_snowpack_steps[:, n_iter] = mbmod.offglac_wide_snowpack
                         output_offglac_runoff_steps[:, n_iter] = mbmod.offglac_wide_runoff
-
-                        if output_glac_bin_icethickness_annual is None:
+                        # binned ouputs
+                        if args.option_dynamics == 'OGGM':
+                            # grab binned outputs from oggm flowline diagnostics
+                            # note, transpose restructures (time, dis_along_flowline) -> (dis_along_flowline, time)
+                            output_glac_bin_area_annual_sim = ds[0].area_m2.values.T[:, :, np.newaxis]
+                            output_glac_bin_icethickness_annual_sim = ds[0].thickness_m.values.T[:, :, np.newaxis]
+                            output_glac_bin_mass_annual_sim = (
+                                ds[0].volume_m3.values.T[:, :, np.newaxis] * pygem_prms['constants']['density_ice']
+                            )
+                            output_glac_bin_flux_divergence_annual_sim = -ds[0].flux_divergence_myr.values.T[
+                                :, :, np.newaxis
+                            ]
+                        else:
                             output_glac_bin_area_annual_sim = mbmod.glac_bin_area_annual[:, :, np.newaxis]
                             output_glac_bin_mass_annual_sim = (
                                 mbmod.glac_bin_area_annual
@@ -1242,6 +1253,9 @@ def run(list_packed_vars):
                             output_glac_bin_icethickness_annual_sim = (mbmod.glac_bin_icethickness_annual)[
                                 :, :, np.newaxis
                             ]
+                            output_glac_bin_flux_divergence_annual_sim = (mbmod.glac_bin_flux_divergence_annual)[
+                                :, :, np.newaxis
+                            ]
                             # Update the latest thickness and volume
                             if ev_model is not None:
                                 fl_dx_meter = getattr(ev_model.fls[0], 'dx_meter', None)
@@ -1263,9 +1277,20 @@ def run(list_packed_vars):
                                 output_glac_bin_mass_annual_sim[:, -1, 0] = (
                                     glacier_vol_t0 * pygem_prms['constants']['density_ice']
                                 )
+                                # flux divergence
+                                output_glac_bin_flux_divergence_annual_sim[:, -1, 0] = (
+                                    mbmod.glac_bin_massbalclim_annual[:, -1]
+                                    * (
+                                        pygem_prms['constants']['density_water']
+                                        / pygem_prms['constants']['density_ice']
+                                    )
+                                ) - np.diff(output_glac_bin_icethickness_annual_sim[:, -2:, 0], axis=-1).flatten()
+                        # append individual simulation outputs together
+                        if output_glac_bin_icethickness_annual is None:
                             output_glac_bin_area_annual = output_glac_bin_area_annual_sim
                             output_glac_bin_mass_annual = output_glac_bin_mass_annual_sim
                             output_glac_bin_icethickness_annual = output_glac_bin_icethickness_annual_sim
+                            output_glac_bin_flux_divergence_annual = output_glac_bin_flux_divergence_annual_sim
                             output_glac_bin_massbalclim_annual_sim = np.zeros(mbmod.glac_bin_icethickness_annual.shape)
                             output_glac_bin_massbalclim_annual_sim[:, :-1] = mbmod.glac_bin_massbalclim_annual
                             output_glac_bin_massbalclim_annual = output_glac_bin_massbalclim_annual_sim[
@@ -1288,36 +1313,6 @@ def run(list_packed_vars):
                             output_glac_bin_melt_steps = output_glac_bin_melt_steps_sim[:, :, np.newaxis]
 
                         else:
-                            # Update the latest thickness and volume
-                            output_glac_bin_area_annual_sim = mbmod.glac_bin_area_annual[:, :, np.newaxis]
-                            output_glac_bin_mass_annual_sim = (
-                                mbmod.glac_bin_area_annual
-                                * mbmod.glac_bin_icethickness_annual
-                                * pygem_prms['constants']['density_ice']
-                            )[:, :, np.newaxis]
-                            output_glac_bin_icethickness_annual_sim = (mbmod.glac_bin_icethickness_annual)[
-                                :, :, np.newaxis
-                            ]
-                            if ev_model is not None:
-                                fl_dx_meter = getattr(ev_model.fls[0], 'dx_meter', None)
-                                fl_widths_m = getattr(ev_model.fls[0], 'widths_m', None)
-                                fl_section = getattr(ev_model.fls[0], 'section', None)
-                            else:
-                                fl_dx_meter = getattr(nfls[0], 'dx_meter', None)
-                                fl_widths_m = getattr(nfls[0], 'widths_m', None)
-                                fl_section = getattr(nfls[0], 'section', None)
-                            if fl_section is not None and fl_widths_m is not None:
-                                # thickness
-                                icethickness_t0 = np.zeros(fl_section.shape)
-                                icethickness_t0[fl_widths_m > 0] = (
-                                    fl_section[fl_widths_m > 0] / fl_widths_m[fl_widths_m > 0]
-                                )
-                                output_glac_bin_icethickness_annual_sim[:, -1, 0] = icethickness_t0
-                                # mass
-                                glacier_vol_t0 = fl_widths_m * fl_dx_meter * icethickness_t0
-                                output_glac_bin_mass_annual_sim[:, -1, 0] = (
-                                    glacier_vol_t0 * pygem_prms['constants']['density_ice']
-                                )
                             output_glac_bin_area_annual = np.append(
                                 output_glac_bin_area_annual,
                                 output_glac_bin_area_annual_sim,
@@ -1333,6 +1328,11 @@ def run(list_packed_vars):
                                 output_glac_bin_icethickness_annual_sim,
                                 axis=2,
                             )
+                            output_glac_bin_flux_divergence_annual = np.append(
+                                output_glac_bin_flux_divergence_annual,
+                                output_glac_bin_flux_divergence_annual_sim,
+                                axis=2,
+                            )
                             output_glac_bin_massbalclim_annual_sim = np.zeros(mbmod.glac_bin_icethickness_annual.shape)
                             output_glac_bin_massbalclim_annual_sim[:, :-1] = mbmod.glac_bin_massbalclim_annual
                             output_glac_bin_massbalclim_annual = np.append(
@@ -1620,6 +1620,9 @@ def run(list_packed_vars):
                                 output_ds_binned_stats['bin_thick_annual'].values[0, :, :] = (
                                     output_glac_bin_icethickness_annual[:, :, n_iter]
                                 )
+                                output_ds_binned_stats['bin_flux_divergence_annual'].values[0, :, :] = (
+                                    output_glac_bin_flux_divergence_annual[:, :, n_iter]
+                                )
                                 output_ds_binned_stats['bin_massbalclim_annual'].values[0, :, :] = (
                                     output_glac_bin_massbalclim_annual[:, :, n_iter]
                                 )
@@ -1679,6 +1682,9 @@ def run(list_packed_vars):
                         output_ds_binned_stats['bin_thick_annual'].values = np.median(
                             output_glac_bin_icethickness_annual, axis=2
                         )[np.newaxis, :, :]
+                        output_ds_binned_stats['bin_flux_divergence_annual'].values = np.median(
+                            output_glac_bin_flux_divergence_annual, axis=2
+                        )[np.newaxis, :, :]
                         output_ds_binned_stats['bin_massbalclim_annual'].values = np.median(
                             output_glac_bin_massbalclim_annual, axis=2
                         )[np.newaxis, :, :]
@@ -1702,6 +1708,9 @@ def run(list_packed_vars):
                             output_ds_binned_stats['bin_thick_annual_mad'].values = median_abs_deviation(
                                 output_glac_bin_icethickness_annual, axis=2
                             )[np.newaxis, :, :]
+                            output_ds_binned_stats['bin_flux_divergence_annual_mad'].values = median_abs_deviation(
+                                output_glac_bin_flux_divergence_annual, axis=2
+                            )[np.newaxis, :, :]
                             output_ds_binned_stats['bin_massbalclim_annual_mad'].values = median_abs_deviation(
                                 output_glac_bin_massbalclim_annual, axis=2
                             )[np.newaxis, :, :]
diff --git a/pygem/massbalance.py b/pygem/massbalance.py
index 20425eb7..e9dd08a7 100644
--- a/pygem/massbalance.py
+++ b/pygem/massbalance.py
@@ -141,6 +141,7 @@ def __init__(
         self.glac_bin_area_annual = np.zeros((nbins, self.nyears + 1))
         self.glac_bin_icethickness_annual = np.zeros((nbins, self.nyears + 1))  # Needed for MassRedistributionCurves
         self.glac_bin_width_annual = np.zeros((nbins, self.nyears + 1))  # Needed for MassRedistributionCurves
+        self.glac_bin_flux_divergence_annual = np.full_like(self.glac_bin_icethickness_annual, np.nan)
         self.offglac_bin_prec = np.zeros((nbins, self.nsteps))
         self.offglac_bin_melt = np.zeros((nbins, self.nsteps))
         self.offglac_bin_refreeze = np.zeros((nbins, self.nsteps))
@@ -763,6 +764,12 @@ def get_annual_mb(
             )
             # Record binned glacier area
             self.glac_bin_area_annual[:, year_idx] = glacier_area_t0
+            # estimate annual flux divergence
+            if year_idx > 0:
+                self.glac_bin_flux_divergence_annual[:, year_idx] = (
+                    self.glac_bin_massbalclim_annual[:, year_idx]
+                    * (pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'])
+                ) - np.diff(self.glac_bin_icethickness_annual[:, year_idx - 1 : year_idx + 1], axis=-1).flatten()
             # Store glacier-wide results
             self._convert_glacwide_results(
                 year_idx,
diff --git a/pygem/output.py b/pygem/output.py
index f1cfa7f8..e806ea50 100644
--- a/pygem/output.py
+++ b/pygem/output.py
@@ -814,6 +814,19 @@ def _update_dicts(self):
             'temporal_resolution': 'annual',
             'comment': 'binned ice thickness at start of the year',
         }
+        self.output_coords_dict['bin_flux_divergence_annual'] = collections.OrderedDict(
+            [
+                ('glac', self.glac_values),
+                ('bin', self.bin_values),
+                ('year', self.year_values),
+            ]
+        )
+        self.output_attrs_dict['bin_flux_divergence_annual'] = {
+            'long_name': 'binned flux divergence, in meters ice',
+            'units': 'm',
+            'temporal_resolution': 'annual',
+            'comment': 'binned flux divergence over the course of the year. note, index i represents the flux divergence over the entire model year[i]',
+        }
         self.output_coords_dict['bin_massbalclim_annual'] = collections.OrderedDict(
             [
                 ('glac', self.glac_values),
@@ -913,6 +926,19 @@ def _update_dicts(self):
                 'temporal_resolution': 'annual',
                 'comment': 'thickness of ice at start of the year',
             }
+            self.output_coords_dict['bin_flux_divergence_annual_mad'] = collections.OrderedDict(
+                [
+                    ('glac', self.glac_values),
+                    ('bin', self.bin_values),
+                    ('year', self.year_values),
+                ]
+            )
+            self.output_attrs_dict['bin_flux_divergence_annual_mad'] = {
+                'long_name': 'binned annual flux divergence median absolute deviation, in meters ice',
+                'units': 'm',
+                'temporal_resolution': 'annual',
+                'comment': 'binned flux divergence over the course of the year. note, index i represents the flux divergence over the entire model year[i]',
+            }
             self.output_coords_dict['bin_massbalclim_annual_mad'] = collections.OrderedDict(
                 [
                     ('glac', self.glac_values),

From 1d4a072eb57eeed2179bc1155881e5057761b55a Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Tue, 11 Nov 2025 10:04:47 -0500
Subject: [PATCH 2/3] Compute annual flux divergence from mb and dhdt

---
 pygem/bin/run/run_calibration.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index eef33cd4..20bae6f3 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -332,9 +332,17 @@ def calc_thick_change_1d(gdir, fls, mbmod, ds):
     """
     years_subannual = np.array([d.year for d in gdir.dates_table['date']])
     yrs = np.unique(years_subannual)
-    nyrs = len(yrs)
     # grab components of interest
     bin_thick_annual = ds[0].thickness_m.values.T  # glacier thickness [m ice], (nbins, nyears)
+    bin_massbalclim_ice_annual = ds[0].climatic_mb_myr.values.T[
+        :, 1:
+    ]  # annual climatic mass balance [m ice] (nbins, nyears - 1)
+    bin_delta_thick_annual = ds[0].dhdt_myr.values.T[
+        :, 1:
+    ]  # annual change in ice thickness [m ice] (nbins, nyears - 1)
+    bin_flux_divergence_annual = (
+        bin_massbalclim_ice_annual - bin_delta_thick_annual
+    )  # annual flux divergence [m ice] (nbins, nyears - 1)
 
     #  --- Step 1: convert mass balance from m w.e. to m ice  ---
     bin_massbalclim = mbmod.glac_bin_massbalclim  # climatic mass balance [m w.e.] per step
@@ -346,7 +354,6 @@ def calc_thick_change_1d(gdir, fls, mbmod, ds):
     # --- Step 2: expand flux divergence to subannual steps ---
     # assume flux divergence is constant throughout the year
     # (divide annual by the number of steps in the binned climatic mass balance to get subannual flux divergence)
-    bin_flux_divergence_annual = -ds[0].flux_divergence_myr.values.T[:, 1:]
     bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice)
     for i, year in enumerate(yrs):
         idx = np.where(years_subannual == year)[0]

From dcc0f7075f3582f23180595710f1c6feb5c313ac Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Tue, 11 Nov 2025 10:36:36 -0500
Subject: [PATCH 3/3] Compute annual flux divergence from mb and dhdt

---
 .../postproc_binned_subannual_thick.py        | 46 ++++++++++---------
 pygem/bin/run/run_simulation.py               | 29 ------------
 pygem/massbalance.py                          |  7 ---
 pygem/output.py                               | 26 -----------
 4 files changed, 25 insertions(+), 83 deletions(-)

diff --git a/pygem/bin/postproc/postproc_binned_subannual_thick.py b/pygem/bin/postproc/postproc_binned_subannual_thick.py
index e917776a..b7c37821 100644
--- a/pygem/bin/postproc/postproc_binned_subannual_thick.py
+++ b/pygem/bin/postproc/postproc_binned_subannual_thick.py
@@ -70,7 +70,6 @@ def get_binned_subannual(
     bin_massbalclim,
     bin_mass_annual,
     bin_thick_annual,
-    bin_flux_divergence_annual,
     dates_subannual,
     dates_annual,
     debug=False,
@@ -85,8 +84,10 @@ def get_binned_subannual(
     thus, subannual thickness and mass is determined assuming
     the flux divergence is constant throughout the year.
 
-    assumes flux divergence is constant throughout the year, dividing annual by the number of steps
-    in the binned climatic mass balance to get subannual flux divergence).
+    annual flux divergence is first calculated by combining the annual binned change in ice
+    thickness and the annual binned mass balance. then, assume flux divergence is constant
+    throughout the year (divide annual by the number of steps in the binned climatic mass
+    balance to get subannual flux divergence).
 
     subannual binned flux divergence can then be combined with
     subannual binned climatic mass balance to get subannual binned change in ice thickness and mass.
@@ -103,9 +104,6 @@ def get_binned_subannual(
     bin_thick_annual : ndarray
        annual binned glacier thickness [m ice]
         shape : [#glac, #elevbins, #years]
-    bin_flux_divergence_annual : ndarray
-       annual binned flux divergence [m ice]
-        shape : [#glac, #elevbins, #years]
     dates_subannual : array-like of datetime-like
         dates associated with `bin_massbalclim` (subannual)
     dates_annual : array-like of datetime-like
@@ -125,35 +123,42 @@ def get_binned_subannual(
         shape : [#glac, #elevbins, #steps]
     """
 
-    n_glac, n_bins, n_steps = bin_massbalclim.shape
-    years_annual = np.array([d.year for d in dates_annual])
     years_subannual = np.array([d.year for d in dates_subannual])
     yrs = np.unique(years_subannual)
-    nyrs = len(yrs)
-    assert nyrs > 1, 'Need at least two annual steps for flux divergence estimation'
+    assert len(yrs) > 1, 'Need at least two annual steps for flux divergence estimation'
 
     #  --- Step 1: convert mass balance from m w.e. to m ice  ---
     rho_w = pygem_prms['constants']['density_water']
     rho_i = pygem_prms['constants']['density_ice']
     bin_massbalclim_ice = bin_massbalclim * (rho_w / rho_i)
 
-    # --- Step 2: expand flux divergence to subannual steps ---
-    # assume flux divergence is constant throughout the year
-    # (divide annual by the number of steps in the binned climatic mass balance to get subannual flux divergence)
-    bin_flux_divergence_annual = bin_flux_divergence_annual[:, :, 1:]
+    # --- Step 2: compute annual thickness change ---
+    bin_delta_thick_annual = np.diff(bin_thick_annual, axis=-1)  # [m ice]
+
+    # --- Step 3: compute annual cumulative mass balance ---
+    # Initialize annual cumulative mass balance (exclude last year for flux calculation)
+    bin_massbalclim_annual = np.zeros_like(bin_delta_thick_annual)
+    for i, year in enumerate(yrs):
+        idx = np.where(years_subannual == year)[0]
+        bin_massbalclim_annual[:, :, i] = bin_massbalclim_ice[:, :, idx].sum(axis=-1)
+
+    # --- Step 4: compute annual flux divergence ---
+    bin_flux_divergence_annual = bin_massbalclim_annual - bin_delta_thick_annual  # [m ice]
+
+    # --- Step 5: expand flux divergence to subannual steps ---
     bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice)
     for i, year in enumerate(yrs):
         idx = np.where(years_subannual == year)[0]
         bin_flux_divergence_subannual[:, :, idx] = bin_flux_divergence_annual[:, :, i][:, :, np.newaxis] / len(idx)
 
-    # --- Step 3: compute subannual thickness change ---
+    # --- Step 6: compute subannual thickness change ---
     bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual
 
-    # --- Step 4: calculate subannual thickness ---
-    running_bin_delta_thick_subannual = np.nancumsum(bin_delta_thick_subannual, axis=-1)
+    # --- Step 7: calculate subannual thickness ---
+    running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1)
     bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, :, 0][:, :, np.newaxis]
 
-    # --- Step 5: compute glacier volume and area on subannual timestep ---
+    # --- Step 8: compute glacier volume and area on subannual timestep ---
     bin_volume_annual = bin_mass_annual / rho_i  # annual volume [m^3] per bin
     bin_area_annual = np.divide(
         bin_volume_annual[:, :, 1:],  # exclude first year to match flux_div_annual
@@ -162,7 +167,7 @@ def get_binned_subannual(
         where=bin_thick_annual[:, :, 1:] > 0,
     )
 
-    # --- Step 6 : compute subannual glacier mass ---
+    # --- Step 9 : compute subannual glacier mass ---
     # First expand area to subannual steps
     bin_area_subannual = np.full(bin_massbalclim_ice.shape, np.nan)
     for i, year in enumerate(yrs):
@@ -172,7 +177,7 @@ def get_binned_subannual(
     # multiply by ice density to get subannual mass
     bin_mass_subannual = bin_thick_subannual * rho_i * bin_area_subannual
 
-    # --- Step 7: debug check ---
+    # --- Step 10: debug check ---
     if debug:
         for i, year in enumerate(yrs):
             # get last subannual index of that year
@@ -299,7 +304,6 @@ def run(simpath, debug):
             bin_massbalclim=binned_ds.bin_massbalclim.values,
             bin_mass_annual=binned_ds.bin_mass_annual.values,
             bin_thick_annual=binned_ds.bin_thick_annual.values,
-            bin_flux_divergence_annual=binned_ds.bin_flux_divergence_annual.values,
             dates_subannual=dates_subannual,
             dates_annual=dates_annual,
             debug=debug,
diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py
index 22ca9935..fafcaf26 100755
--- a/pygem/bin/run/run_simulation.py
+++ b/pygem/bin/run/run_simulation.py
@@ -1240,9 +1240,6 @@ def run(list_packed_vars):
                             output_glac_bin_mass_annual_sim = (
                                 ds[0].volume_m3.values.T[:, :, np.newaxis] * pygem_prms['constants']['density_ice']
                             )
-                            output_glac_bin_flux_divergence_annual_sim = -ds[0].flux_divergence_myr.values.T[
-                                :, :, np.newaxis
-                            ]
                         else:
                             output_glac_bin_area_annual_sim = mbmod.glac_bin_area_annual[:, :, np.newaxis]
                             output_glac_bin_mass_annual_sim = (
@@ -1253,9 +1250,6 @@ def run(list_packed_vars):
                             output_glac_bin_icethickness_annual_sim = (mbmod.glac_bin_icethickness_annual)[
                                 :, :, np.newaxis
                             ]
-                            output_glac_bin_flux_divergence_annual_sim = (mbmod.glac_bin_flux_divergence_annual)[
-                                :, :, np.newaxis
-                            ]
                             # Update the latest thickness and volume
                             if ev_model is not None:
                                 fl_dx_meter = getattr(ev_model.fls[0], 'dx_meter', None)
@@ -1277,20 +1271,11 @@ def run(list_packed_vars):
                                 output_glac_bin_mass_annual_sim[:, -1, 0] = (
                                     glacier_vol_t0 * pygem_prms['constants']['density_ice']
                                 )
-                                # flux divergence
-                                output_glac_bin_flux_divergence_annual_sim[:, -1, 0] = (
-                                    mbmod.glac_bin_massbalclim_annual[:, -1]
-                                    * (
-                                        pygem_prms['constants']['density_water']
-                                        / pygem_prms['constants']['density_ice']
-                                    )
-                                ) - np.diff(output_glac_bin_icethickness_annual_sim[:, -2:, 0], axis=-1).flatten()
                         # append individual simulation outputs together
                         if output_glac_bin_icethickness_annual is None:
                             output_glac_bin_area_annual = output_glac_bin_area_annual_sim
                             output_glac_bin_mass_annual = output_glac_bin_mass_annual_sim
                             output_glac_bin_icethickness_annual = output_glac_bin_icethickness_annual_sim
-                            output_glac_bin_flux_divergence_annual = output_glac_bin_flux_divergence_annual_sim
                             output_glac_bin_massbalclim_annual_sim = np.zeros(mbmod.glac_bin_icethickness_annual.shape)
                             output_glac_bin_massbalclim_annual_sim[:, :-1] = mbmod.glac_bin_massbalclim_annual
                             output_glac_bin_massbalclim_annual = output_glac_bin_massbalclim_annual_sim[
@@ -1328,11 +1313,6 @@ def run(list_packed_vars):
                                 output_glac_bin_icethickness_annual_sim,
                                 axis=2,
                             )
-                            output_glac_bin_flux_divergence_annual = np.append(
-                                output_glac_bin_flux_divergence_annual,
-                                output_glac_bin_flux_divergence_annual_sim,
-                                axis=2,
-                            )
                             output_glac_bin_massbalclim_annual_sim = np.zeros(mbmod.glac_bin_icethickness_annual.shape)
                             output_glac_bin_massbalclim_annual_sim[:, :-1] = mbmod.glac_bin_massbalclim_annual
                             output_glac_bin_massbalclim_annual = np.append(
@@ -1620,9 +1600,6 @@ def run(list_packed_vars):
                                 output_ds_binned_stats['bin_thick_annual'].values[0, :, :] = (
                                     output_glac_bin_icethickness_annual[:, :, n_iter]
                                 )
-                                output_ds_binned_stats['bin_flux_divergence_annual'].values[0, :, :] = (
-                                    output_glac_bin_flux_divergence_annual[:, :, n_iter]
-                                )
                                 output_ds_binned_stats['bin_massbalclim_annual'].values[0, :, :] = (
                                     output_glac_bin_massbalclim_annual[:, :, n_iter]
                                 )
@@ -1682,9 +1659,6 @@ def run(list_packed_vars):
                         output_ds_binned_stats['bin_thick_annual'].values = np.median(
                             output_glac_bin_icethickness_annual, axis=2
                         )[np.newaxis, :, :]
-                        output_ds_binned_stats['bin_flux_divergence_annual'].values = np.median(
-                            output_glac_bin_flux_divergence_annual, axis=2
-                        )[np.newaxis, :, :]
                         output_ds_binned_stats['bin_massbalclim_annual'].values = np.median(
                             output_glac_bin_massbalclim_annual, axis=2
                         )[np.newaxis, :, :]
@@ -1708,9 +1682,6 @@ def run(list_packed_vars):
                             output_ds_binned_stats['bin_thick_annual_mad'].values = median_abs_deviation(
                                 output_glac_bin_icethickness_annual, axis=2
                             )[np.newaxis, :, :]
-                            output_ds_binned_stats['bin_flux_divergence_annual_mad'].values = median_abs_deviation(
-                                output_glac_bin_flux_divergence_annual, axis=2
-                            )[np.newaxis, :, :]
                             output_ds_binned_stats['bin_massbalclim_annual_mad'].values = median_abs_deviation(
                                 output_glac_bin_massbalclim_annual, axis=2
                             )[np.newaxis, :, :]
diff --git a/pygem/massbalance.py b/pygem/massbalance.py
index e9dd08a7..20425eb7 100644
--- a/pygem/massbalance.py
+++ b/pygem/massbalance.py
@@ -141,7 +141,6 @@ def __init__(
         self.glac_bin_area_annual = np.zeros((nbins, self.nyears + 1))
         self.glac_bin_icethickness_annual = np.zeros((nbins, self.nyears + 1))  # Needed for MassRedistributionCurves
         self.glac_bin_width_annual = np.zeros((nbins, self.nyears + 1))  # Needed for MassRedistributionCurves
-        self.glac_bin_flux_divergence_annual = np.full_like(self.glac_bin_icethickness_annual, np.nan)
         self.offglac_bin_prec = np.zeros((nbins, self.nsteps))
         self.offglac_bin_melt = np.zeros((nbins, self.nsteps))
         self.offglac_bin_refreeze = np.zeros((nbins, self.nsteps))
@@ -764,12 +763,6 @@ def get_annual_mb(
             )
             # Record binned glacier area
             self.glac_bin_area_annual[:, year_idx] = glacier_area_t0
-            # estimate annual flux divergence
-            if year_idx > 0:
-                self.glac_bin_flux_divergence_annual[:, year_idx] = (
-                    self.glac_bin_massbalclim_annual[:, year_idx]
-                    * (pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'])
-                ) - np.diff(self.glac_bin_icethickness_annual[:, year_idx - 1 : year_idx + 1], axis=-1).flatten()
             # Store glacier-wide results
             self._convert_glacwide_results(
                 year_idx,
diff --git a/pygem/output.py b/pygem/output.py
index e806ea50..f1cfa7f8 100644
--- a/pygem/output.py
+++ b/pygem/output.py
@@ -814,19 +814,6 @@ def _update_dicts(self):
             'temporal_resolution': 'annual',
             'comment': 'binned ice thickness at start of the year',
         }
-        self.output_coords_dict['bin_flux_divergence_annual'] = collections.OrderedDict(
-            [
-                ('glac', self.glac_values),
-                ('bin', self.bin_values),
-                ('year', self.year_values),
-            ]
-        )
-        self.output_attrs_dict['bin_flux_divergence_annual'] = {
-            'long_name': 'binned flux divergence, in meters ice',
-            'units': 'm',
-            'temporal_resolution': 'annual',
-            'comment': 'binned flux divergence over the course of the year. note, index i represents the flux divergence over the entire model year[i]',
-        }
         self.output_coords_dict['bin_massbalclim_annual'] = collections.OrderedDict(
             [
                 ('glac', self.glac_values),
@@ -926,19 +913,6 @@ def _update_dicts(self):
                 'temporal_resolution': 'annual',
                 'comment': 'thickness of ice at start of the year',
             }
-            self.output_coords_dict['bin_flux_divergence_annual_mad'] = collections.OrderedDict(
-                [
-                    ('glac', self.glac_values),
-                    ('bin', self.bin_values),
-                    ('year', self.year_values),
-                ]
-            )
-            self.output_attrs_dict['bin_flux_divergence_annual_mad'] = {
-                'long_name': 'binned annual flux divergence median absolute deviation, in meters ice',
-                'units': 'm',
-                'temporal_resolution': 'annual',
-                'comment': 'binned flux divergence over the course of the year. note, index i represents the flux divergence over the entire model year[i]',
-            }
             self.output_coords_dict['bin_massbalclim_annual_mad'] = collections.OrderedDict(
                 [
                     ('glac', self.glac_values),