Merge 2fdf052 into 9ddfb51

docs/whats-new.rst

@@ -32,6 +32,15 @@ Enhancements
   By `Anouk Vlug <https://github.com/anoukvlug>`_ and
   `Fabien Maussion <https://github.com/fmaussion>`_
 
+Enhancements
+~~~~~~~~~~~~
+
+- Added three new flowline diagnostic variables: thickness change in one year
+  (``dhdt``), forcing climatic mass balance (``climatic_mb``) and flux divergence
+  (``flux_divergence``). All variables are in units meter of ice per year
+  (:pull:`1595`).
+  By `Patrick Schmitt <https://github.com/pat-schmitt>`_
+
 
 v1.6.0 (March 10, 2023)
 -----------------------

oggm/core/flowline.py

@@ -1214,6 +1214,44 @@ def run_until_and_store(self, y1,
                     desc = 'Part of the series which are spinup'
                     ds['is_fixed_geometry_spinup'].attrs['description'] = desc
                     ds['is_fixed_geometry_spinup'].attrs['unit'] = '-'
+                if 'flux_divergence' in ovars_fl:
+                    # We need dhdt and climatic_mb to calculate divergence
+                    if 'dhdt' not in ovars_fl:
+                        ovars_fl.append('dhdt')
+                    if 'climatic_mb' not in ovars_fl:
+                        ovars_fl.append('climatic_mb')
+
+                    ds['flux_divergence_myr'] = (('time', 'dis_along_flowline'),
+                                                 width * np.NaN)
+                    desc = 'Ice-flux divergence'
+                    ds['flux_divergence_myr'].attrs['description'] = desc
+                    ds['flux_divergence_myr'].attrs['unit'] = 'm yr-1'
+                if 'climatic_mb' in ovars_fl:
+                    # We need dhdt to define where to calculate climatic_mb
+                    if 'dhdt' not in ovars_fl:
+                        ovars_fl.append('dhdt')
+
+                    ds['climatic_mb_myr'] = (('time', 'dis_along_flowline'),
+                                             width * np.NaN)
+                    desc = 'Climatic mass balance forcing'
+                    ds['climatic_mb_myr'].attrs['description'] = desc
+                    ds['climatic_mb_myr'].attrs['unit'] = 'm yr-1'
+
+                    # needed for the calculation of mb at start of period
+                    surface_h_previous = {}
+                    for fl_id, fl in enumerate(self.fls):
+                        surface_h_previous[fl_id] = fl.surface_h
+                if 'dhdt' in ovars_fl:
+                    ds['dhdt_myr'] = (('time', 'dis_along_flowline'),
+                                      width * np.NaN)
+                    desc = 'Thickness change'
+                    ds['dhdt_myr'].attrs['description'] = desc
+                    ds['dhdt_myr'].attrs['unit'] = 'm yr-1'
+
+                    # needed for the calculation of dhdt
+                    thickness_previous_dhdt = {}
+                    for fl_id, fl in enumerate(self.fls):
+                        thickness_previous_dhdt[fl_id] = fl.thick
 
         # First deal with spinup (we compute volume change only)
         if do_fixed_spinup:
@@ -1269,6 +1307,47 @@ def run_until_and_store(self, y1,
                             var = self.u_stag[fl_id]
                             val = (var[1:fl.nx + 1] + var[:fl.nx]) / 2 * self._surf_vel_fac
                             ds['ice_velocity_myr'].data[j, :] = val * cfg.SEC_IN_YEAR
+                        if 'dhdt' in ovars_fl and (yr > self.y0):
+                            # dhdt can only be computed after one year
+                            val = fl.thick - thickness_previous_dhdt[fl_id]
+                            ds['dhdt_myr'].data[j, :] = val
+                            thickness_previous_dhdt[fl_id] = fl.thick
+                        if 'climatic_mb' in ovars_fl and (yr > self.y0):
+                            # yr - 1 to use the mb which lead to the current
+                            # state, also using previous surface height
+                            val = self.get_mb(surface_h_previous[fl_id],
+                                              self.yr - 1,
+                                              fl_id=fl_id)
+                            # only save climatic mb where dhdt is non zero,
+                            # isclose for avoiding numeric represention artefacts
+                            dhdt_zero = np.isclose(ds['dhdt_myr'].data[j, :],
+                                                   0.)
+                            ds['climatic_mb_myr'].data[j, :] = np.where(
+                                dhdt_zero,
+                                0.,
+                                val * cfg.SEC_IN_YEAR)
+                            surface_h_previous[fl_id] = fl.surface_h
+                        if 'flux_divergence' in ovars_fl and (yr > self.y0):
+                            # calculated after the formula dhdt = mb + flux_div
+                            val = ds['dhdt_myr'].data[j, :] -\
+                                  ds['climatic_mb_myr'].data[j, :]
+                            # special treatment for retreating: If the glacier
+                            # terminus is getting ice free it means the
+                            # climatic mass balance is more negative than the
+                            # available ice to melt. In this case we can not
+                            # calculate the flux divergence with this method
+                            # exactly, we just can say it was smaller to
+                            # compensate the very negative climatic mb.
+                            # To avoid large spikes at the terminus in the
+                            # calculated flux divergence we smooth this values
+                            # with an (arbitrary) factor of 0.1.
+                            has_become_ice_free = np.logical_and(
+                                np.isclose(fl.thick, 0.),
+                                ds['dhdt_myr'].data[j, :] < 0.
+                            )
+                            fac = np.where(has_become_ice_free, 0.1, 1.)
+                            ds['flux_divergence_myr'].data[j, :] = val * fac
+
                 # j is the yearly index in case we have monthly output
                 # we have to count it ourselves
                 j += 1

oggm/params.cfg

@@ -420,4 +420,4 @@ store_fl_diagnostics = False
 # What variables would you like to store in the flowline diagnostics files
 # Note: area and volume are mandatory
 # You need to keep all variables in one line unfortunately
-store_fl_diagnostic_variables =  area, thickness, volume, volume_bsl, volume_bwl, calving_bucket, ice_velocity
+store_fl_diagnostic_variables =  area, thickness, volume, volume_bsl, volume_bwl, calving_bucket, ice_velocity, dhdt, climatic_mb, flux_divergence

oggm/tests/test_models.py

@@ -1612,6 +1612,11 @@ def test_run(self, class_case_dir):
         vol_ref = []
         a_ref = []
         l_ref = []
+        h_previous_timestep = model.fls[0].thick
+        dhdt_ref = []
+        surface_h_previous_timestep = model.fls[0].surface_h
+        climatic_mb_ref = []
+        flux_divergence_ref = []
         vol_diag = []
         a_diag = []
         l_diag = []
@@ -1626,6 +1631,29 @@ def test_run(self, class_case_dir):
                 vol_ref.append(model.volume_m3)
                 a_ref.append(model.area_m2)
                 l_ref.append(model.length_m)
+                if yr > 0:
+                    dhdt_ref.append(model.fls[0].thick -
+                                    h_previous_timestep)
+                    h_previous_timestep = model.fls[0].thick
+                    # for mb use previous surface height and previous year, only
+                    # save climatic mb where dhdt is non zero
+                    climatic_mb_ref.append(
+                        np.where(np.isclose(dhdt_ref[-1], 0.),
+                                 0.,
+                                 model.get_mb(surface_h_previous_timestep,
+                                              model.yr - 1,
+                                              fl_id=0) *
+                                 SEC_IN_YEAR)  # converted to m yr-1
+                    )
+                    surface_h_previous_timestep = model.fls[0].surface_h
+                    # smooth flux divergence where glacier is getting ice free
+                    has_become_ice_free = np.logical_and(
+                                np.isclose(model.fls[0].thick, 0.),
+                                dhdt_ref[-1] < 0.
+                            )
+                    flux_divergence_ref.append(
+                        (dhdt_ref[-1] - climatic_mb_ref[-1]) *
+                        np.where(has_become_ice_free, 0.1, 1.))
                 if int(yr) == 500:
                     secfortest = model.fls[0].section
                     hfortest = model.fls[0].thick
@@ -1648,6 +1676,13 @@ def test_run(self, class_case_dir):
         np.testing.assert_allclose(ds_fl.area_m2.sum(dim='dis_along_flowline'),
                                    a_ref)
 
+        np.testing.assert_allclose(dhdt_ref,
+                                   ds_fl.dhdt_myr[1:])  # first step is nan
+        np.testing.assert_allclose(climatic_mb_ref,
+                                   ds_fl.climatic_mb_myr[1:])
+        np.testing.assert_allclose(flux_divergence_ref,
+                                   ds_fl.flux_divergence_myr[1:])
+
         vel = ds_fl.ice_velocity_myr.isel(time=-1)
         assert 20 < vel.max() < 40