diff --git a/docs/whats-new.rst b/docs/whats-new.rst
index c044e831f..edbd631e7 100644
--- a/docs/whats-new.rst
+++ b/docs/whats-new.rst
@@ -55,6 +55,8 @@ Enhancements
   By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 - Added support for Millan et al 2022 velocity and thickness in the shop (:pull:`1443`).
   By `Fabien Maussion <https://github.com/fmaussion>`_
+- Added support for Hugonnet et al 2021 dhdt in the shop (:pull:`1529`).
+  By `Fabien Maussion <https://github.com/fmaussion>`_
 - Added trapezoidal downstream line (:pull:`1491`). Can be selected with
   ``cfg.PARAMS['downstream_line_shape']``, with the options ``'parabol'`` (default)
   or ``'trapezoidal'`` before calling ``init_present_time_glacier(gdir)``.
diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py
index 6383b7264..510cbcf6b 100644
--- a/oggm/cli/prepro_levels.py
+++ b/oggm/cli/prepro_levels.py
@@ -81,8 +81,10 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
                       match_geodetic_mb_per_glacier=False,
                       evolution_model='fl_sia',
                       centerlines_only=False, override_params=None,
-                      add_consensus=False, start_level=None,
-                      start_base_url=None, max_level=5, ref_tstars_base_url='',
+                      add_consensus_thickness=False, add_millan_thickness=False,
+                      add_millan_velocity=False, add_hugonnet_dhdt=False,
+                      start_level=None, start_base_url=None, max_level=5,
+                      ref_tstars_base_url='',
                       logging_level='WORKFLOW', disable_dl_verify=False,
                       dynamic_spinup=False, err_dmdtda_scaling_factor=1,
                       dynamic_spinup_start_year=1979,
@@ -139,9 +141,18 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
     evolution_model : str
         which geometry evolution model to use: `fl_sia` (default),
         or `massredis` (mass redistribution curve).
-    add_consensus : bool
+    add_consensus_thickness : bool
         adds (reprojects) the consensus estimates thickness to the glacier
         directories. With elev_bands=True, the data will also be binned.
+    add_millan_thickness : bool
+        adds (reprojects) the millan thickness to the glacier
+        directories. With elev_bands=True, the data will also be binned.
+    add_millan_velocity : bool
+        adds (reprojects) the millan velocity to the glacier
+        directories. With elev_bands=True, the data will also be binned.
+    add_hugonnet_dhdt : bool
+        adds (reprojects) the hugonnet dhdt maps to the glacier
+        directories. With elev_bands=True, the data will also be binned.
     start_level : int
         the pre-processed level to start from (default is to start from
         scratch). If set, you'll need to indicate start_base_url as well.
@@ -430,17 +441,33 @@ def _time_log():
         # Centerlines OGGM
         workflow.execute_entity_task(tasks.glacier_masks, gdirs_cent)
 
-        if add_consensus:
+        bin_variables = []
+        if add_consensus_thickness:
             from oggm.shop.bedtopo import add_consensus_thickness
             workflow.execute_entity_task(add_consensus_thickness, gdirs_band)
             workflow.execute_entity_task(add_consensus_thickness, gdirs_cent)
-
-            # Elev bands with var data
-            vn = 'consensus_ice_thickness'
+            bin_variables.append('consensus_ice_thickness')
+        if add_millan_thickness:
+            from oggm.shop.millan22 import thickness_to_gdir
+            workflow.execute_entity_task(thickness_to_gdir, gdirs_band)
+            workflow.execute_entity_task(thickness_to_gdir, gdirs_cent)
+            bin_variables.append('millan_ice_thickness')
+        if add_millan_velocity:
+            from oggm.shop.millan22 import velocity_to_gdir
+            workflow.execute_entity_task(velocity_to_gdir, gdirs_band)
+            workflow.execute_entity_task(velocity_to_gdir, gdirs_cent)
+            bin_variables.append('millan_v')
+        if add_hugonnet_dhdt:
+            from oggm.shop.hugonnet_maps import hugonnet_to_gdir
+            workflow.execute_entity_task(hugonnet_to_gdir, gdirs_band)
+            workflow.execute_entity_task(hugonnet_to_gdir, gdirs_cent)
+            bin_variables.append('hugonnet_dhdt')
+
+        if bin_variables and gdirs_band:
             workflow.execute_entity_task(tasks.elevation_band_flowline,
-                                         gdirs_band, bin_variables=vn)
+                                         gdirs_band, bin_variables=bin_variables)
             workflow.execute_entity_task(tasks.fixed_dx_elevation_band_flowline,
-                                         gdirs_band, bin_variables=vn)
+                                         gdirs_band, bin_variables=bin_variables)
         else:
             # HH2015 method without it
             task_list = [
@@ -480,6 +507,19 @@ def _time_log():
         opath = os.path.join(sum_dir, 'glacier_statistics_{}.csv'.format(rgi_reg))
         utils.compile_glacier_statistics(gdirs, path=opath)
 
+        if add_millan_thickness or add_millan_velocity:
+            from oggm.shop.millan22 import compile_millan_statistics
+            opath = os.path.join(sum_dir, 'millan_statistics_{}.csv'.format(rgi_reg))
+            compile_millan_statistics(gdirs, path=opath)
+        if add_hugonnet_dhdt:
+            from oggm.shop.hugonnet_maps import compile_hugonnet_statistics
+            opath = os.path.join(sum_dir, 'hugonnet_statistics_{}.csv'.format(rgi_reg))
+            compile_hugonnet_statistics(gdirs, path=opath)
+        if add_consensus_thickness:
+            from oggm.shop.bedtopo import compile_consensus_statistics
+            opath = os.path.join(sum_dir, 'consensus_statistics_{}.csv'.format(rgi_reg))
+            compile_consensus_statistics(gdirs, path=opath)
+
         # And for level 2: shapes
         if len(gdirs_cent) > 0:
             opath = os.path.join(sum_dir, f'centerlines_{rgi_reg}.shp')
@@ -814,9 +854,24 @@ def parse_args(args):
                              'compatible with level 1 folders, after which '
                              'the processing will stop. The default is to use '
                              'the default OGGM DEM.')
-    parser.add_argument('--add-consensus', nargs='?', const=True, default=False,
-                        help='adds (reprojects) the consensus estimates '
-                             'thickness to the glacier directories. '
+    parser.add_argument('--add-consensus-thickness', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the consensus thickness '
+                             'estimates to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-millan-thickness', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the millan thickness '
+                             'estimates to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-millan-velocity', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the millan velocity '
+                             'estimates to the glacier directories. '
+                             'With --elev-bands, the data will also be '
+                             'binned.')
+    parser.add_argument('--add-hugonnet-dhdt', nargs='?', const=True, default=False,
+                        help='adds (reprojects) the millan dhdt '
+                             'maps to the glacier directories. '
                              'With --elev-bands, the data will also be '
                              'binned.')
     parser.add_argument('--demo', nargs='?', const=True, default=False,
@@ -899,7 +954,10 @@ def parse_args(args):
                 centerlines_only=args.centerlines_only,
                 match_regional_geodetic_mb=args.match_regional_geodetic_mb,
                 match_geodetic_mb_per_glacier=args.match_geodetic_mb_per_glacier,
-                add_consensus=args.add_consensus,
+                add_consensus_thickness=args.add_consensus_thickness,
+                add_millan_thickness=args.add_millan_thickness,
+                add_millan_velocity=args.add_millan_velocity,
+                add_hugonnet_dhdt=args.add_hugonnet_dhdt,
                 disable_dl_verify=args.disable_dl_verify,
                 ref_tstars_base_url=args.ref_tstars_base_url,
                 evolution_model=args.evolution_model,
diff --git a/oggm/shop/bedtopo.py b/oggm/shop/bedtopo.py
index 376322cb7..b9c0301d4 100644
--- a/oggm/shop/bedtopo.py
+++ b/oggm/shop/bedtopo.py
@@ -1,13 +1,16 @@
 import logging
+import os
 
 import numpy as np
+import pandas as pd
+import xarray as xr
 
 try:
     import salem
 except ImportError:
     pass
 
-from oggm import utils
+from oggm import utils, cfg
 
 # Module logger
 log = logging.getLogger(__name__)
@@ -68,3 +71,62 @@ def add_consensus_thickness(gdir, base_url=None):
         v.long_name = ln
         v.base_url = base_url
         v[:] = thick
+
+
+@utils.entity_task(log)
+def consensus_statistics(gdir):
+    """Gather statistics about the consensus data interpolated to this glacier.
+    """
+
+    d = dict()
+
+    # Easy stats - this should always be possible
+    d['rgi_id'] = gdir.rgi_id
+    d['rgi_region'] = gdir.rgi_region
+    d['rgi_subregion'] = gdir.rgi_subregion
+    d['rgi_area_km2'] = gdir.rgi_area_km2
+    d['consensus_vol_km3'] = 0
+    d['consensus_area_km2'] = 0
+    d['consensus_perc_cov'] = 0
+
+    try:
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+            thick = ds['consensus_ice_thickness'].where(ds['glacier_mask'], np.NaN).load()
+            d['consensus_vol_km3'] = float(thick.sum() * gdir.grid.dx ** 2 * 1e-9)
+            d['consensus_area_km2'] = float((~thick.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6)
+            d['consensus_perc_cov'] = float(d['consensus_area_km2'] / gdir.rgi_area_km2)
+    except (FileNotFoundError, AttributeError, KeyError):
+        pass
+
+    return d
+
+
+@utils.global_task(log)
+def compile_consensus_statistics(gdirs, filesuffix='', path=True):
+    """Gather as much statistics as possible about a list of glaciers.
+
+    Parameters
+    ----------
+    gdirs : list of :py:class:`oggm.GlacierDirectory` objects
+        the glacier directories to process
+    filesuffix : str
+        add suffix to output file
+    path : str, bool
+        Set to "True" in order  to store the info in the working directory
+        Set to a path to store the file to your chosen location
+    """
+    from oggm.workflow import execute_entity_task
+
+    out_df = execute_entity_task(consensus_statistics, gdirs)
+
+    out = pd.DataFrame(out_df).set_index('rgi_id')
+
+    if path:
+        if path is True:
+            out.to_csv(os.path.join(cfg.PATHS['working_dir'],
+                                    ('consensus_statistics' +
+                                     filesuffix + '.csv')))
+        else:
+            out.to_csv(path)
+
+    return out
diff --git a/oggm/shop/hugonnet_maps.py b/oggm/shop/hugonnet_maps.py
new file mode 100644
index 000000000..26c831763
--- /dev/null
+++ b/oggm/shop/hugonnet_maps.py
@@ -0,0 +1,211 @@
+import logging
+from packaging.version import Version
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+import os
+
+try:
+    import rasterio
+    from rasterio.warp import reproject, Resampling
+    from rasterio import MemoryFile
+    try:
+        # rasterio V > 1.0
+        from rasterio.merge import merge as merge_tool
+    except ImportError:
+        from rasterio.tools.merge import merge as merge_tool
+except ImportError:
+    pass
+
+try:
+    import salem
+except ImportError:
+    pass
+
+try:
+    import geopandas as gpd
+except ImportError:
+    pass
+
+from oggm import utils, cfg
+
+# Module logger
+log = logging.getLogger(__name__)
+
+default_base_url = 'https://cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb_maps/'
+
+_lookup_csv = None
+
+
+def _get_lookup_csv():
+    global _lookup_csv
+    if _lookup_csv is None:
+        fname = default_base_url + 'hugonnet_dhdt_lookup_csv_20230129.csv'
+        _lookup_csv = pd.read_csv(utils.file_downloader(fname), index_col=0)
+    return _lookup_csv
+
+
+@utils.entity_task(log, writes=['gridded_data'])
+def hugonnet_to_gdir(gdir, add_error=False):
+    """Add the Hugonnet 21 dhdt maps to this glacier directory.
+
+    Parameters
+    ----------
+    gdir : :py:class:`oggm.GlacierDirectory`
+        the glacier directory to process
+    add_error : bool
+        add the error data or not
+    """
+
+    if add_error:
+        raise NotImplementedError('Not yet')
+
+    # Find out which file(s) we need
+    df = _get_lookup_csv()
+    lon_ex, lat_ex = gdir.extent_ll
+
+    # adding small buffer for unlikely case where one lon/lat_ex == xx.0
+    lons = np.arange(np.floor(lon_ex[0] - 1e-9), np.ceil(lon_ex[1] + 1e-9))
+    lats = np.arange(np.floor(lat_ex[0] - 1e-9), np.ceil(lat_ex[1] + 1e-9))
+
+    flist = []
+    for lat in lats:
+        # north or south?
+        ns = 'S' if lat < 0 else 'N'
+        for lon in lons:
+            # east or west?
+            ew = 'W' if lon < 0 else 'E'
+            ll_str = f'{ns}{abs(lat):02.0f}{ew}{abs(lon):03.0f}'
+            try:
+                filename = df.loc[(df['file_id'] == ll_str)]['dhdt'].iloc[0]
+            except IndexError:
+                # We can maybe be on the edge (unlikely but hey
+                pass
+            file_local = utils.file_downloader(default_base_url + filename)
+            if file_local is not None:
+                flist.append(file_local)
+
+    # A glacier area can cover more than one tile:
+    if len(flist) == 1:
+        dem_dss = [rasterio.open(flist[0])]  # if one tile, just open it
+        dem_data = rasterio.band(dem_dss[0], 1)
+        if Version(rasterio.__version__) >= Version('1.0'):
+            src_transform = dem_dss[0].transform
+        else:
+            src_transform = dem_dss[0].affine
+        nodata = dem_dss[0].meta.get('nodata', None)
+    else:
+        dem_dss = [rasterio.open(s) for s in flist]  # list of rasters
+        nodata = dem_dss[0].meta.get('nodata', None)
+        dem_data, src_transform = merge_tool(dem_dss, nodata=nodata)  # merge
+
+    # Set up profile for writing output
+    with rasterio.open(gdir.get_filepath('dem')) as dem_ds:
+        dst_array = dem_ds.read().astype(np.float32)
+        dst_array[:] = np.NaN
+        profile = dem_ds.profile
+        transform = dem_ds.transform
+        dst_crs = dem_ds.crs
+
+    # Set up profile for writing output
+    profile.update({
+        'nodata': np.NaN,
+    })
+
+    resampling = Resampling.bilinear
+
+    with MemoryFile() as dest:
+        reproject(
+            # Source parameters
+            source=dem_data,
+            src_crs=dem_dss[0].crs,
+            src_transform=src_transform,
+            src_nodata=nodata,
+            # Destination parameters
+            destination=dst_array,
+            dst_transform=transform,
+            dst_crs=dst_crs,
+            dst_nodata=np.NaN,
+            # Configuration
+            resampling=resampling)
+        dest.write(dst_array)
+
+    for dem_ds in dem_dss:
+        dem_ds.close()
+
+    # Write
+    with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc:
+
+        vn = 'hugonnet_dhdt'
+        if vn in nc.variables:
+            v = nc.variables[vn]
+        else:
+            v = nc.createVariable(vn, 'f4', ('y', 'x', ), zlib=True)
+        v.units = 'm'
+        ln = 'dhdt (2000-2020) from Hugonnet et al. 2021'
+        v.long_name = ln
+        data_str = ' '.join(flist) if len(flist) > 1 else flist[0]
+        v.data_source = data_str
+        v[:] = np.squeeze(dst_array)
+
+
+@utils.entity_task(log)
+def hugonnet_statistics(gdir):
+    """Gather statistics about the Hugonnet data interpolated to this glacier.
+    """
+
+    d = dict()
+
+    # Easy stats - this should always be possible
+    d['rgi_id'] = gdir.rgi_id
+    d['rgi_region'] = gdir.rgi_region
+    d['rgi_subregion'] = gdir.rgi_subregion
+    d['rgi_area_km2'] = gdir.rgi_area_km2
+    d['hugonnet_area_km2'] = 0
+    d['hugonnet_perc_cov'] = 0
+    d['hugonnet_avg_dhdt'] = np.NaN
+
+    try:
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+            dhdt = ds['hugonnet_dhdt'].where(ds['glacier_mask'], np.NaN).load()
+            d['hugonnet_area_km2'] = float((~dhdt.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6)
+            d['hugonnet_perc_cov'] = float(d['hugonnet_area_km2'] / gdir.rgi_area_km2)
+            d['hugonnet_avg_dhdt'] = np.nanmean(dhdt.data)
+    except (FileNotFoundError, AttributeError, KeyError):
+        pass
+
+    return d
+
+
+@utils.global_task(log)
+def compile_hugonnet_statistics(gdirs, filesuffix='', path=True):
+    """Gather as much statistics as possible about a list of glaciers.
+
+    It can be used to do result diagnostics and other stuffs.
+
+    Parameters
+    ----------
+    gdirs : list of :py:class:`oggm.GlacierDirectory` objects
+        the glacier directories to process
+    filesuffix : str
+        add suffix to output file
+    path : str, bool
+        Set to "True" in order  to store the info in the working directory
+        Set to a path to store the file to your chosen location
+    """
+    from oggm.workflow import execute_entity_task
+
+    out_df = execute_entity_task(hugonnet_statistics, gdirs)
+
+    out = pd.DataFrame(out_df).set_index('rgi_id')
+
+    if path:
+        if path is True:
+            out.to_csv(os.path.join(cfg.PATHS['working_dir'],
+                                    ('hugonnet_statistics' +
+                                     filesuffix + '.csv')))
+        else:
+            out.to_csv(path)
+
+    return out
diff --git a/oggm/tests/test_shop.py b/oggm/tests/test_shop.py
index 7a2e6c1d8..d18fd627e 100644
--- a/oggm/tests/test_shop.py
+++ b/oggm/tests/test_shop.py
@@ -1,6 +1,7 @@
 import os
 import warnings
 
+import matplotlib.pyplot as plt
 import pytest
 salem = pytest.importorskip('salem')
 gpd = pytest.importorskip('geopandas')
@@ -12,7 +13,7 @@
 import pandas as pd
 from oggm import utils
 from oggm.utils import get_demo_file
-from oggm.shop import its_live, rgitopo, bedtopo, millan22
+from oggm.shop import its_live, rgitopo, bedtopo, millan22, hugonnet_maps
 from oggm.core import gis, centerlines, massbalance
 from oggm import cfg, tasks, workflow
 
@@ -160,6 +161,43 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
         assert df['millan_vol_km3'] > 174
 
 
+class Test_HugonnetMaps:
+
+    @pytest.mark.slow
+    def test_dhdt_to_glacier(self, class_case_dir, monkeypatch):
+
+        # Init
+        cfg.initialize()
+        cfg.PATHS['working_dir'] = class_case_dir
+        cfg.PARAMS['use_intersects'] = False
+        cfg.PATHS['dem_file'] = get_demo_file('dem_Columbia.tif')
+        cfg.PARAMS['border'] = 10
+
+        entity = gpd.read_file(get_demo_file('RGI60-01.10689.shp')).iloc[0]
+        gdir = oggm.GlacierDirectory(entity)
+        tasks.define_glacier_region(gdir)
+        tasks.glacier_masks(gdir)
+
+        # use our files
+        base_url = 'https://cluster.klima.uni-bremen.de/~oggm/test_files/geodetic_ref_mb_maps/'
+        monkeypatch.setattr(hugonnet_maps, 'default_base_url', base_url)
+
+        hugonnet_maps.hugonnet_to_gdir(gdir)
+
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+            mask = ds.glacier_mask.data.astype(bool)
+            dhdt = ds.hugonnet_dhdt.where(mask)
+
+        # Simply some coverage and sanity checks
+        assert np.isfinite(dhdt).sum() / mask.sum() > 0.97
+        assert np.nanmax(dhdt) > 0
+        assert np.nanmean(dhdt) < -3
+
+        df = hugonnet_maps.compile_hugonnet_statistics([gdir]).iloc[0]
+        assert df['hugonnet_perc_cov'] > 0.97
+        assert df['hugonnet_avg_dhdt'] < -3
+
+
 class Test_rgitopo:
 
     def test_from_dem(self, class_case_dir, monkeypatch):
@@ -768,3 +806,7 @@ def test_add_consensus(self, class_case_dir, monkeypatch):
         # Check vol
         my_vol = (fdf[vn] * fdf['area_m2']).sum()
         np.testing.assert_allclose(my_vol, ref_vol)
+
+        df = bedtopo.compile_consensus_statistics([gdir]).iloc[0]
+        np.testing.assert_allclose(my_vol*1e-9, df['consensus_vol_km3'], rtol=1e-2)
+        np.testing.assert_allclose(1, df['consensus_perc_cov'], rtol=0.05)
diff --git a/oggm/tests/test_utils.py b/oggm/tests/test_utils.py
index 507e30a3a..b4795ea40 100644
--- a/oggm/tests/test_utils.py
+++ b/oggm/tests/test_utils.py
@@ -864,6 +864,7 @@ def test_parse_args(self):
         assert kwargs['border'] == 160
         assert not kwargs['dynamic_spinup']
         assert kwargs['dynamic_spinup_start_year'] == 1979
+        assert not kwargs['add_consensus_thickness']
 
         kwargs = prepro_levels.parse_args(['--rgi-reg', '1',
                                            '--map-border', '160',
@@ -931,6 +932,10 @@ def test_parse_args(self):
                                            '--output', 'local/out',
                                            '--working-dir', 'local/work',
                                            '--dem-source', 'ALL',
+                                           '--add-consensus-thickness',
+                                           '--add-millan-thickness',
+                                           '--add-millan-velocity',
+                                           '--add-hugonnet-dhdt',
                                            ])
 
         assert 'working_dir' in kwargs
@@ -945,6 +950,10 @@ def test_parse_args(self):
         assert not kwargs['elev_bands']
         assert not kwargs['match_regional_geodetic_mb']
         assert not kwargs['centerlines_only']
+        assert kwargs['add_consensus_thickness']
+        assert kwargs['add_millan_thickness']
+        assert kwargs['add_millan_velocity']
+        assert kwargs['add_hugonnet_dhdt']
 
         kwargs = prepro_levels.parse_args(['--rgi-reg', '1',
                                            '--map-border', '160',