diff --git a/src/lisflood/Lisflood_dynamic.py b/src/lisflood/Lisflood_dynamic.py
index e5ffe47d..77d15892 100644
--- a/src/lisflood/Lisflood_dynamic.py
+++ b/src/lisflood/Lisflood_dynamic.py
@@ -246,6 +246,10 @@ def splitlanduse(array1, array2=None, array3=None):
         # Calculate water level
         self.waterlevel_module.dynamic()
 
+        # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        # Calculate water storage
+        self.waterstorage_module.dynamic()        
+
         # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
         # ************************************************************
diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py
index 346c1dff..697aef74 100644
--- a/src/lisflood/Lisflood_initial.py
+++ b/src/lisflood/Lisflood_initial.py
@@ -59,6 +59,7 @@
 from .hydrological_modules.opensealed import opensealed
 from .hydrological_modules.waterbalance import waterbalance
 from .hydrological_modules.waterlevel import waterlevel
+from .hydrological_modules.waterstorage import waterstorage
 from .hydrological_modules.structures import structures
 
 from .global_modules.output import outputTssMap
@@ -150,6 +151,7 @@ def __init__(self):
         self.opensealed_module = opensealed(self)
         self.waterbalance_module = waterbalance(self)
         self.waterlevel_module = waterlevel(self)
+        self.waterstorage_module = waterstorage(self)
         self.structures_module = structures(self)
 
         self.prescribed_vegetation = self.epic_settings.prescribed_vegetation
@@ -204,6 +206,7 @@ def __init__(self):
         self.reservoir_module.initial()
         self.lakes_module.initial()
         self.polder_module.initial()
+        self.waterstorage_module.initial()
 
         self.transmission_module.initial()
 
diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py
index 3091665f..86139337 100644
--- a/src/lisflood/global_modules/default_options.py
+++ b/src/lisflood/global_modules/default_options.py
@@ -1597,3 +1597,12 @@
                    'wateruseRegion': False,
                    'writeNetcdf': False,
                    'writeNetcdfStack': False}
+
+default_options['reportedmaps'].update({'TWSMaps'     : ReportedMap(name='TWSMaps',    output_var='twsstor',       unit='m',     end=[], steps=['repTWSMaps'],                      all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'LakeSMaps'   : ReportedMap(name='LakeSMaps',  output_var='lakestor',      unit='m3',    end=[], steps=['repStorageMaps'],                  all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'RiverSMaps'  : ReportedMap(name='RiverSMaps', output_var='riverstor',     unit='m3',    end=[], steps=['repStorageMaps','repFlowMomMaps'], all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'SoilSMaps'   : ReportedMap(name='SoilSMaps',  output_var='soilstor',      unit='m',     end=[], steps=['repStorageMaps'],                  all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'GWSMaps'     : ReportedMap(name='GWSMaps',    output_var='gwstor',        unit='m',     end=[], steps=['repStorageMaps'],                  all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'SnowSMaps'   : ReportedMap(name='SnowSMaps',  output_var='snowstor',      unit='m',     end=[], steps=['repStorageMaps'],                  all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'CumSMaps '   : ReportedMap(name='CumSMaps',   output_var='cumstor',       unit='m',     end=[], steps=['repStorageMaps'],                  all=[], restrictoption=[], monthly=False, yearly=False)})
+default_options['reportedmaps'].update({'FlowMomMaps' : ReportedMap(name='FlowMomMaps', output_var='FlowMomentum', unit='kgm/s', end=[], steps=['repFlowMomMaps'],                  all=[], restrictoption=[], monthly=False, yearly=False)})
diff --git a/src/lisflood/hydrological_modules/evapowater.py b/src/lisflood/hydrological_modules/evapowater.py
index a7cc18db..ef67fff1 100755
--- a/src/lisflood/hydrological_modules/evapowater.py
+++ b/src/lisflood/hydrological_modules/evapowater.py
@@ -60,6 +60,7 @@ def initial(self):
         maskinfo = MaskInfo.instance()
         if option['openwaterevapo']:
             LakeMask = loadmap('LakeMask', pcr=True)
+            LakeMask = ifthenelse((LakeMask<0)&(LakeMask>-9999), 0, LakeMask)
             # lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresKinematic, 5)
             lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresChan, 5)
             LddEva = lddrepair(lmask)
diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py
index f4d0dcba..c0b7624f 100644
--- a/src/lisflood/hydrological_modules/routing.py
+++ b/src/lisflood/hydrological_modules/routing.py
@@ -931,6 +931,11 @@ def dynamic(self, NoRoutingExecuted):
             # if flow is slow, Traveltime=DtSec then TravelDistance=PixelLength
             # maximum set to 30km/day for 5km cell, is at DtSec/Traveltime=6, is at Traveltime<DtSec/6
 
+            ####### added by Robert Dill, 2025/03/13
+            # calculate momentum map flow velocity * rivermass [kgm/s]
+            self.var.FlowMomentum = self.var.TotalCrossSectionArea * self.var.ChanLength * 1000 * self.var.FlowVelocity
+            #######
+
 
     def KinRouting(self,SideflowChan):
         """Kinematic routing based on a 4-point implicit finite-difference numerical solution of the kinematic wave equations.
diff --git a/src/lisflood/hydrological_modules/waterstorage.py b/src/lisflood/hydrological_modules/waterstorage.py
new file mode 100644
index 00000000..ac6a5c26
--- /dev/null
+++ b/src/lisflood/hydrological_modules/waterstorage.py
@@ -0,0 +1,296 @@
+"""
+
+Copyright 2019 European Union
+
+Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
+
+You may not use this work except in compliance with the Licence.
+You may obtain a copy of the Licence at:
+
+https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
+
+Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the Licence for the specific language governing permissions and limitations under the Licence.
+
+"""
+from __future__ import absolute_import, print_function
+
+import warnings
+warnings.formatwarning = lambda msg, args, *kwargs: f'{msg}\n'
+
+import numpy as np
+from ..global_modules.add1 import loadmap
+from ..global_modules.settings import LisSettings, MaskInfo
+from ..global_modules.errors import LisfloodError, LisfloodWarning
+from . import HydroModule
+
+
+class waterstorage(HydroModule):
+    """
+    # ************************************************************
+    # ***** WATER STORAGE    *************************************
+    # ************************************************************
+    # Sum up water storage in individual compartements 
+    """
+    input_files_keys = {
+        'repTWSMaps': ['LakeMask'],
+        'repStorageMaps': ['LakeMask']
+    }
+    module_name = 'WaterStorage'
+
+    def __init__(self, waterstorage_variable):
+        self.var = waterstorage_variable
+
+    # --------------------------------------------------------------------------
+    # --------------------------------------------------------------------------
+
+    def initial(self):
+        """ initial part of the water storage module
+        """
+
+        # ************************************************************
+        # ***** WATER STORAGE INIT
+        # ************************************************************
+        settings = LisSettings.instance()
+        option = settings.options
+        
+        # load water map and separate into maps of lake and reservoir distribution
+        if (not(option['InitLisflood'])) and (option['repStorageMaps'] or option['repTWSMaps']):
+
+            LakeMask = loadmap('LakeMask')
+            # find ranges of IDs for lakes [LakeID_min+1:LakeID_max] and reservoirs [ReservoirID_min+1:ReservoirID_max]
+            # LakeMask   land       = 0
+            #            water      = 1
+            #                         or for lakes:      ID = LakeID_min(default:1000)+lakeID
+            #                                reservoirs: ID = ReservoirID_min(default:5000) + reservoirID
+            #                         in case LakeMask==0 (no water defined before insertion of lake-/reservoir-IDs)
+            #                         lake and reservoirs IDs are inserted with negative values
+            #
+            # the gap between 1 and lakeIDs and between lakeIds and reservoirIDs has to be greater than <gap_min>
+            # gaps within lakeIDs or reservoirIDs have to be smaller than <gap_min>
+            gap_min = 499
+            LakeMask[LakeMask<=-9999] = np.nan
+            LakeMask = np.abs(LakeMask)
+            id_all = np.unique(LakeMask[~np.isnan(LakeMask)]).astype(int)
+            id_max = id_all[-1]
+            lowerBounds = (id_all+1)[:-1]
+            upperBounds = (id_all-1)[1:]
+            mask = lowerBounds<=upperBounds-(max(1,gap_min)-1)
+
+            n_gaps = np.count_nonzero(mask)
+            if n_gaps == 0:
+                # only land/water mask, no lakes or reservoir IDs
+                LakeID_min = id_max
+                LakeID_max = LakeID_min
+                ReservoirID_min = id_max
+                ReservoirID_max = ReservoirID_min
+            else:
+                # beginning of lake IDs
+                LakeID_min = upperBounds[mask][0]
+                if n_gaps == 1:
+                    # only lake IDs
+                    LakeID_max = id_max+1
+                    # no reservoir IDs
+                    ReservoirID_min = id_max
+                    ReservoirID_max = ReservoirID_min
+                else:
+                    # end of lake IDs
+                    LakeID_max = lowerBounds[mask][1]
+                    # beginning of reservoir IDs
+                    ReservoirID_min = upperBounds[mask][1]
+                    # end of reservoir IDs
+                    ReservoirID_max = id_max+1
+
+            # extract lake distribution map
+            self.var.LakeDistribution = np.where((LakeMask > LakeID_min) & (LakeMask < LakeID_max), LakeMask-LakeID_min, np.nan)
+            # extract reservoir distribution map
+            self.var.ReservoirDistribution = np.where((LakeMask > ReservoirID_min) & (LakeMask < ReservoirID_max), LakeMask-ReservoirID_min, np.nan)
+
+            # check number of ID with number of sites
+            if n_gaps < 1:
+                msg = "{} LakeMask map (containing no lake or reservoir IDs) not compatible for TWS calculation"
+                raise LisfloodError(msg)
+            if option['simulateLakes']:
+                number_of_lakeIDs = len(np.unique(self.var.LakeDistribution[~np.isnan(self.var.LakeDistribution)]).astype(int))
+                if self.var.LakeSitesCC.size != number_of_lakeIDs:
+                    warnings.warn(LisfloodWarning('Number of lake IDs in map LakeMask ('+str(number_of_lakeIDs)+') not equal number of lake sites defined in map LakeSites ('+str(self.var.LakeSitesCC.size)+').'))
+            if option['simulateReservoirs']:
+                number_of_reservoirIDs = len(np.unique(self.var.ReservoirDistribution[~np.isnan(self.var.ReservoirDistribution)]).astype(int))
+                if self.var.ReservoirSitesCC.size != number_of_reservoirIDs:
+                    warnings.warn(LisfloodWarning('Number of reservoir IDs in map LakeMask ('+str(number_of_reservoirIDs)+') not equal number of reservoir sites defined in map ReservoirSites ('+str(self.var.ReservoirSitesCC.size)+').'))
+
+
+# --------------------------------------------------------------------------
+# -------------------------------------------------------------------------
+
+    def dynamic(self):
+        """ dynamic part of the water storage module
+        """
+        settings = LisSettings.instance()
+        option = settings.options
+        maskinfo = MaskInfo.instance()
+            
+        
+        if (not(option['InitLisflood'])) and (option['repStorageMaps'] or option['repTWSMaps']):
+
+            # ************************************************************
+            # ***** WATER STORAGE 
+            # ************************************************************
+            
+            # river water storage [m3]
+            # in routing.py: *** initial river storage is stored in
+            #                    self.var.ChanM3 = self.var.TotalCrossSectionArea * self.var.ChanLength
+            #                    self.var.ChanIniM3 = self.var.ChanM3.copy()
+            #                    self.var.ChanM3Kin = self.var.ChanIniM3.copy().astype(float)
+            #                *** split routing
+            #                    self.var.Chan2M3Kin = self.var.CrossSection2Area * self.var.ChanLength + self.var.Chan2M3Start
+            #                    self.var.ChanM3Kin = self.var.ChanM3 - self.var.Chan2M3Kin + self.var.Chan2M3Start
+            #                *** Volume in main channel at end of computation step
+            #                    self.var.ChanM3Kin = self.var.ChanLength * self.var.ChannelAlpha * self.var.ChanQKin**self.var.Beta
+            #                *** floodplain routing
+            #                    self.var.Chan2M3Kin = self.var.ChanLength * self.var.ChannelAlpha2 * self.var.Chan2QKin ** self.var.Beta
+            #                    self.var.CrossSection2Area = (self.var.Chan2M3Kin - self.var.Chan2M3Start) * self.var.InvChanLength  
+            #                    TotalCrossSectionArea = np.maximum(self.var.ChanM3Kin*self.var.InvChanLength,0.01)
+            # in Lisflood_dynamic.py:
+            #                *** add main channel and floodplains
+            #                    self.ChanM3 = self.ChanM3Kin + selfChan2M3Kin - self.Chan2M3Start
+            #                    self.TotalCrossSectionArea = self.ChanM3 * self.InvChanLength            
+            #tws_riverM3 = self.var.ChanM3.copy() ?
+            tws_riverM3 = self.var.TotalCrossSectionArea * self.var.ChanLength
+            
+            # [m3] -> [m]
+            tws_riverM = tws_riverM3 / self.var.PixelArea
+            
+            # overlandflow water storage [m3]
+            # in surface_routing.py: self.var.M3all = self.var.OFM3Direct + self.var.OFM3Other + self.var.OFM3Forest
+            tws_oflowM3 = self.var.OFM3Direct + self.var.OFM3Forest + self.var.OFM3Other
+            
+            # [m3] -> [m]
+            tws_oflowM = tws_oflowM3 / self.var.PixelArea
+            
+            # lake water storage [m3]
+            # in lakes.py: LakeSitesC = loadmap('LakeSites')
+            #              LakeSitesC[LakeSitesC < 1] = 0
+            #              LakeSitesC[self.var.IsChannel == 0] = 0
+            #              self.var.LakeSitesC2 = LakeSitesC 
+            #
+            #              LakeArea = pcraster.lookupscalar(str(binding['TabLakeArea']), LakeSitePcr)
+            #              LakeAreaC = compressArray(LakeArea)
+            #              self.var.LakeAreaCC = np.compress(LakeSitesC > 0, LakeAreaC)
+            #
+            #              self.var.LakeStorageM3CC = (LakeStorageIndicator - self.var.LakeOutflowCC* 0.5) * self.var.DtRouting
+            #              self.var.LakeStorageM3CC[self.var.LakeStorageM3CC < 0] = 0
+            #              self.var.LakeStorageM3CC[np.isnan(self.var.LakeStorageM3CC)] = 0
+            #              self.var.LakeStorageM3BalanceCC += LakeIn * self.var.DtRouting - QLakeOutM3DtCC
+            #              self.var.LakeLevelCC = self.var.LakeStorageM3CC / self.var.LakeAreaCC
+            #
+            #              self.var.LakeStorageM3Balance = maskinfo.in_zero()
+            #              self.var.LakeStorageM3 = maskinfo.in_zero()
+            #              self.var.LakeLevel = maskinfo.in_zero()
+            #              np.put(self.var.LakeStorageM3Balance, self.var.LakeIndex, self.var.LakeStorageM3BalanceCC)
+            #              np.put(self.var.LakeStorageM3, self.var.LakeIndex, self.var.LakeStorageM3CC)
+            #              np.put(self.var.LakeLevel, self.var.LakeIndex, self.var.LakeLevelCC)
+            tws_lakeM3 = np.zeros(tws_riverM3.shape, dtype=np.float32)
+            if option['simulateLakes']:
+                #tws_lakeM3  = self.var.LakeStorageM3.copy() ? 
+                #tws_lakeM3  = self.var.LakeStorageM3Balance ?
+                LakeArea = maskinfo.in_zero()
+                np.put(LakeArea, self.var.LakeIndex, self.var.LakeAreaCC)
+                tws_lakeM3 = self.var.LakeLevel * LakeArea
+                
+            # [m3] -> [m] distribute lake/river/oflow over lake areas 
+            tws_lakeM = np.zeros(tws_riverM.shape, dtype=np.float32)
+            lake_extent = self.var.LakeDistribution
+            for n in np.unique(lake_extent[~np.isnan(lake_extent)]).astype(int):
+                if n > 0:
+                    lake_mask = np.nonzero(lake_extent == n)
+                    grid_area_lake = np.nansum(self.var.PixelArea[lake_mask])
+                    tws_lakeM[lake_mask] = tws_lakeM3[self.var.LakeSitesC2==n] / grid_area_lake
+                    tws_riverM[lake_mask] = np.nansum(tws_riverM3[lake_mask]) / grid_area_lake
+                    tws_oflowM[lake_mask] = np.nansum(tws_oflowM3[lake_mask]) / grid_area_lake
+            
+            # reservoir water storage [m3]
+            # in routing.py: self.var.IsChannelPcr = boolean(loadmap('Channels', pcr=True))
+            #                self.var.IsChannel = np.bool8(compressArray(self.var.IsChannelPcr))
+            #
+            # in reservoir.py: self.var.ReservoirSitesC = loadmap('ReservoirSites')
+            #                  self.var.ReservoirSitesC[self.var.ReservoirSitesC < 1] = 0
+            #                  self.var.ReservoirSitesC[self.var.IsChannel == 0] = 0
+            #
+            #                  TotalReservoirStorageM3 = lookupscalar(str(binding['TabTotStorage']), ReservoirSitePcr)
+            #                  self.var.TotalReservoirStorageM3C = compressArray(TotalReservoirStorageM3)
+            #                  self.var.TotalReservoirStorageM3C = np.where(np.isnan(self.var.TotalReservoirStorageM3C), 0, self.var.TotalReservoirStorageM3C)
+            #                  self.var.TotalReservoirStorageM3CC = np.compress(self.var.ReservoirSitesC > 0, self.var.TotalReservoirStorageM3C)
+            #
+            #                  self.var.ReservoirStorageM3CC -= QResOutM3DtCC
+            #                  self.var.ReservoirFillCC = self.var.ReservoirStorageM3CC / self.var.TotalReservoirStorageM3CC
+            #                  self.var.ReservoirFillCC[np.isnan(self.var.ReservoirFillCC)] = 0
+            #                  self.var.ReservoirFillCC[self.var.ReservoirFillCC < 0] = 0
+            #
+            #                  self.var.ReservoirStorageM3 = maskinfo.in_zero()
+            #                  self.var.ReservoirFill = maskinfo.in_zero()
+            #                  np.put(self.var.ReservoirStorageM3, self.var.ReservoirIndex, self.var.ReservoirStorageM3CC)
+            #                  np.put(self.var.ReservoirFill, self.var.ReservoirIndex, self.var.ReservoirFillCC)
+            tws_reservoirM3 = np.zeros(tws_riverM3.shape, dtype=np.float32)
+            if option['simulateReservoirs']:  
+                #tws_reservoirM3  = self.var.ReservoirStorageM3.copy() ?
+                TotalReservoirStorage = maskinfo.in_zero()
+                np.put(TotalReservoirStorage, self.var.ReservoirIndex, self.var.TotalReservoirStorageM3CC)
+                tws_reservoirM3 = self.var.ReservoirFill * TotalReservoirStorage
+            
+            # [m3] -> [m] distribute reservoir/river/oflow over reservoir areas 
+            tws_reservoirM = np.zeros(tws_riverM.shape, dtype=np.float32)
+            reservoir_extent = self.var.ReservoirDistribution
+            for n in np.unique(reservoir_extent[~np.isnan(reservoir_extent)]).astype(int):
+                if n > 0:
+                    reservoir_mask = np.nonzero(reservoir_extent == n)
+                    grid_area_reservoir = np.nansum(self.var.PixelArea[reservoir_mask])
+                    tws_reservoirM[reservoir_mask] = tws_reservoirM3[self.var.ReservoirSitesC==n] / grid_area_reservoir
+                    tws_riverM[reservoir_mask] = np.nansum(tws_riverM3[reservoir_mask]) / grid_area_reservoir
+                    tws_oflowM[reservoir_mask] = np.nansum(tws_oflowM3[reservoir_mask]) / grid_area_reservoir
+         
+            # soil water storage [mm] -> [m]
+            tws_soil1M  = ((self.var.Theta1a[0] * self.var.SoilDepth1a[0]) * self.var.OtherFraction     +
+                          (self.var.Theta1a[1] * self.var.SoilDepth1a[1]) * self.var.ForestFraction    +
+                          (self.var.Theta1a[2] * self.var.SoilDepth1a[2]) * self.var.IrrigationFraction
+                         ) / 1000
+
+            tws_soil2M  = ((self.var.Theta1b[0] * self.var.SoilDepth1b[0]) * self.var.OtherFraction     +
+                         (self.var.Theta1b[1] * self.var.SoilDepth1b[1]) * self.var.ForestFraction    +
+                         (self.var.Theta1b[2] * self.var.SoilDepth1b[2]) * self.var.IrrigationFraction
+                        ) / 1000
+
+            tws_soil3M  = ((self.var.Theta2[0] * self.var.SoilDepth2[0]) * self.var.OtherFraction     +
+                          (self.var.Theta2[1] * self.var.SoilDepth2[1]) * self.var.ForestFraction    +
+                          (self.var.Theta2[2] * self.var.SoilDepth2[2]) * self.var.IrrigationFraction
+                        ) / 1000
+
+            tws_soilM = tws_soil1M + tws_soil2M + tws_soil3M
+                                  
+            # groundwater storage [mm] -> [m], (uz,uzf,uzi,lz,fracforest,fracirrigated,fracother)
+            tws_groundwaterM = ( self.var.UZ[0] * self.var.OtherFraction      + 
+                                 self.var.UZ[1] * self.var.ForestFraction     +
+                                 self.var.UZ[2] * self.var.IrrigationFraction +
+                                 self.var.LZ  
+                               ) / 1000
+                               
+            # snow storage [m]
+            tws_snowM = self.var.SnowCover / 1000
+            
+            # cumulative interception storage [mm] and cumulative depression storage [mm] -> [m]
+            tws_cumM = ( self.var.CumInterception[0] * self.var.OtherFraction       +
+                         self.var.CumInterception[1] * self.var.ForestFraction      +
+                         self.var.CumInterception[2] * self.var.IrrigationFraction  +
+                         self.var.CumInterSealed  * self.var.DirectRunoffFraction
+                       ) / 1000
+
+            # total water storage maps
+            self.var.riverstor     = tws_riverM + tws_oflowM
+            self.var.lakestor      = tws_lakeM + tws_reservoirM
+            self.var.soilstor      = tws_soilM
+            self.var.gwstor        = tws_groundwaterM
+            self.var.snowstor      = tws_snowM
+            self.var.cumstor       = tws_cumM
+            self.var.twsstor       = tws_riverM + tws_oflowM + tws_lakeM + tws_reservoirM + tws_soilM + tws_groundwaterM + tws_snowM + tws_cumM
+                      
diff --git a/src/lisfloodSettings_reference.xml b/src/lisfloodSettings_reference.xml
index 297c42d4..8c4d4d19 100644
--- a/src/lisfloodSettings_reference.xml
+++ b/src/lisfloodSettings_reference.xml
@@ -161,6 +161,10 @@ You can use builtin path variables in this template and reference to other paths
       <setoption choice="0" name="repTotalWUse"/>
       <setoption choice="0" name="repWIndex"/>
 
+      <setoption choice="1" name="repTWSMaps"/>
+      <setoption choice="1" name="repStorageMaps"/>
+      <setoption choice="1" name="repFlowMomMaps"/>
+
   </lfoptions>
   <lfuser>
 
@@ -3276,6 +3280,62 @@ Multiplier []  applied to ChanSdXdY
       Reported transpiration maps, weighted sum over the fractions of each pixel [mm]
       </comment>
       </textvar>
+
+      <textvar name="TWSMaps" value="$(PathOut)/TWS">
+      <comment>
+      Reported total water storage storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="LakeSMaps" value="$(PathOut)/Lakestor">
+      <comment>
+      Reported lake and reservoir storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="RiverSMaps" value="$(PathOut)/Riverstor">
+      <comment>
+      Reported river storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="SoilSMaps" value="$(PathOut)/Soilstor">
+      <comment>
+      Reported soil storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="GWSMaps" value="$(PathOut)/GWstor">
+      <comment>
+      Reported groundwater storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="SnowSMaps" value="$(PathOut)/Snowstor">
+      <comment>
+      Reported snow storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="CumSMaps" value="$(PathOut)/Cumstor">
+      <comment>
+      Reported cumulative interception and depression storage storage [m]
+      added by GFZ
+      </comment>
+      </textvar>
+
+      <textvar name="FlowMomMaps" value="$(PathOut)/FlowMomentum">
+      <comment>
+      Reported riverflow momentum for HAM motion term [kgm/s]
+      added by GFZ
+      </comment>
+      </textvar>        
       </group>
 
 
@@ -6166,4 +6226,4 @@ Multiplier []  applied to ChanSdXdY
     ReportSteps=$(ReportSteps);
   </prolog>
 
-</lfsettings>
\ No newline at end of file
+</lfsettings>