Merge pull request #146 from ec-jrc/development

Total water storage reporting

VERSION

@@ -1 +1 @@
-4.3.0
+4.3.1

docs/3_step3_preparing-setting-file/lisfloodSettings_reference.xml

@@ -3378,14 +3378,6 @@ soil layer 2 (lower layer) [V/V] [mm3/mm3]
 </comment>
 </textvar>
 
-<textvar name="ThetaMaps" value="$(PathOut)/th">
-<comment>
-Reported volumetric soil moisture content for  combined
-soil layer 1 and 2 [V/V] [mm3/mm3]
-</comment> 
-</textvar>
-
-
 <textvar name="UZMaps" value="$(PathOut)/uz">
 <comment>
 Reported storage in upper response box [mm]
@@ -3510,7 +3502,7 @@ soil layer 2 (lower layer) [V/V] [mm3/mm3]
 </comment>
 </textvar>
 
-<textvar name="UZForestMaps" value="$(PathOut)/uzF">
+<textvar name="UZForestMaps" value="$(PathOut)/uzf">
 <comment>
 Reported storage in upper response box [mm]
 </comment>
@@ -3783,6 +3775,12 @@ Reported upper zone outflow [mm]
 </comment>
 </textvar>
 
+<textvar name="UZOutflowIrrigationMaps" value="$(PathOut)/quzi">
+<comment>
+Reported upper zone outflow [mm]
+</comment>
+</textvar>
+
 <textvar name="GwPercUZLZForestMaps" value="$(PathOut)/uz2lzF">
 <comment>
 Reported percolation from upper to lower zone [mm]

src/lisflood/__init__.py

@@ -10,7 +10,7 @@
 
 __version__ = version
 __authors__ = "Ad de Roo, Emiliano Gelati, Peter Burek, Johan van der Knijff, Niko Wanders"
-__date__ = "05/01/2024"
+__date__ = "23/01/2024"
 __copyright__ = "Copyright 2019-2024, European Commission - Joint Research Centre"
 __maintainer__ = "Stefania Grimaldi, Cinzia Mazzetti, Carlo Russo, Damien Decremer, Corentin Carton De Wiart, Valerio Lorini, Ad de Roo"
 __status__ = "Operation"

src/lisflood/global_modules/default_options.py

@@ -51,6 +51,7 @@
                    'repLZOutflowMaps': False,
                    'repLeafDrainageMaps': False,
                    'repMBTs': False,
+                   'repTotalWaterStorageMaps': False,
                    'repMeteoUpsGauges': False,
                    'repPFMaps': False,
                    'repPFSites': False,
@@ -222,6 +223,12 @@
                                                       end=['repEndMaps'], steps=[], all=[],
                                                       restrictoption=[], monthly=False,
                                                       yearly=False),
+                       'CumIntSealedMaps': ReportedMap(name='CumIntSealedMaps',
+                                                       output_var='CumInterSealed',
+                                                       unit='mm', end=[], steps=[],
+                                                       all=['repCumInterCeptionMaps'],
+                                                       restrictoption=['nonInit'],
+                                                       monthly=False, yearly=False),
                        'CumIntSealedState': ReportedMap(name='CumIntSealedState',
                                                         output_var='CumInterSealed', unit='mm',
                                                         end=[], steps=['repStateMaps'], all=[],
@@ -232,6 +239,17 @@
                                                          end=['repEndMaps'], steps=[], all=[],
                                                          restrictoption=[], monthly=False,
                                                          yearly=False),
+                       'CumInterceptionMaps': ReportedMap(name='CumInterceptionMaps',
+                                                          output_var='CumInterception[0]', unit='mm',
+                                                          end=[], steps=[],
+                                                          all=['repCumInterCeptionMaps'],
+                                                          restrictoption=['nonInit'], monthly=False,
+                                                          yearly=False),
+                       'CumInterceptionState': ReportedMap(name='CumInterceptionState',
+                                                           output_var='CumInterception[0]', unit='mm',
+                                                           end=[], steps=['repStateMaps'], all=[],
+                                                           restrictoption=['nonInit'], monthly=False,
+                                                           yearly=False),
                        'CumInterceptionForestEnd': ReportedMap(name='CumInterceptionForestEnd',
                                                                output_var='CumInterception[1]',
                                                                unit='mm', end=['repEndMaps'],
@@ -256,21 +274,16 @@
                                                                    steps=[], all=[],
                                                                    restrictoption=[],
                                                                    monthly=False, yearly=False),
+                       'CumInterceptionIrrigationMaps': ReportedMap(name='CumInterceptionIrrigationMaps',
+                                                                    output_var='CumInterception[2]',
+                                                                    unit='mm', end=[], steps=[],
+                                                                    all=['repCumInterCeptionMaps'],
+                                                                    restrictoption=['nonInit'],
+                                                                    monthly=False, yearly=False),
                        'CumInterceptionIrrigationState': ReportedMap(
                            name='CumInterceptionIrrigationState', output_var='CumInterception[2]',
                            unit='mm', end=[], steps=['repStateMaps'], all=[],
                            restrictoption=['nonInit'], monthly=False, yearly=False),
-                       'CumInterceptionMaps': ReportedMap(name='CumInterceptionMaps',
-                                                          output_var='CumInterception[0]', unit='mm',
-                                                          end=[], steps=[],
-                                                          all=['repCumInterCeptionMaps'],
-                                                          restrictoption=['nonInit'], monthly=False,
-                                                          yearly=False),
-                       'CumInterceptionState': ReportedMap(name='CumInterceptionState',
-                                                           output_var='CumInterception[0]', unit='mm',
-                                                           end=[], steps=['repStateMaps'], all=[],
-                                                           restrictoption=['nonInit'], monthly=False,
-                                                           yearly=False),
                        'DSLREnd': ReportedMap(name='DSLREnd', output_var='DSLR[0]', unit='-',
                                               end=['repEndMaps'], steps=[], all=[],
                                               restrictoption=[], monthly=False, yearly=False),
@@ -1099,6 +1112,10 @@
                                                       unit='mm', end=['repEndMaps'], steps=[],
                                                       all=[], restrictoption=[],
                                                       monthly=False, yearly=False),
+                       'UZIrrigationMaps': ReportedMap(name='UZIrrigationMaps', output_var='UZ[2]',
+                                                       unit='mm', end=[], steps=[],
+                                                       all=['repUZMaps'], restrictoption=['nonInit'],
+                                                       monthly=False, yearly=False),
                        'UZIrrigationState': ReportedMap(name='UZIrrigationState', output_var='UZ[2]',
                                                         unit='mm', end=[], steps=['repStateMaps'],
                                                         all=[], restrictoption=['nonInit'],
@@ -1167,6 +1184,12 @@
                                                       unit='m', end=[], steps=[],
                                                       all=[], restrictoption=['nonInit'],
                                                       monthly=False, yearly=False),
+                       'TotalWaterStorageMaps': ReportedMap(name='TotalWaterStorageMaps',
+                                                            output_var='TotalWaterStorageMM',
+                                                            unit='mm', end=[], steps=[],
+                                                            all=['repTotalWaterStorageMaps'],
+                                                            restrictoption=['nonInit'],
+                                                            monthly=False, yearly=False),
                        'WaterSecurityIndex': ReportedMap(name='WaterSecurityIndex',
                                                          output_var='WaterSecurityIndex', unit='-',
                                                          end=[], steps=['repWIndex'], all=[],

src/lisflood/hydrological_modules/waterbalance.py

@@ -108,6 +108,24 @@ def initial(self):
 
             self.var.DischargeM3StructuresIni = np.take(np.bincount(self.var.Catchments, weights=DisStructure), self.var.Catchments)
 
+# --------------------------------------------------------------------------
+# -------------------------------------------------------------------------
+
+    def storage_channel(self, option):
+        ChannelStoredM3 = self.var.ChanM3.copy()
+        if option['simulateLakes']:
+            ChannelStoredM3 += self.var.LakeStorageM3Balance
+        if option['simulateReservoirs']:
+            ChannelStoredM3 += self.var.ReservoirStorageM3
+        if option['simulatePolders']:
+            ChannelStoredM3 += self.var.PolderStorageM3
+        return ChannelStoredM3
+
+    def storage_hillslope(self):
+        ax_veg = self.var.SoilFraction.dims.index("vegetation")
+        Hill1 = self.var.LZ + np.sum(self.var.SoilFraction * (self.var.CumInterception + self.var.W1 + self.var.W2 + self.var.UZ), ax_veg)
+        HillslopeStoredMM = self.var.WaterDepth + self.var.SnowCover + Hill1 + self.var.DirectRunoffFraction * self.var.CumInterSealed
+        return HillslopeStoredMM * self.var.MMtoM3
 
 # --------------------------------------------------------------------------
 # -------------------------------------------------------------------------
@@ -118,6 +136,11 @@ def dynamic(self):
         settings = LisSettings.instance()
         option = settings.options        
         maskinfo = MaskInfo.instance()
+
+        if (not(option['InitLisflood'])) and option['repTotalWaterStorageMaps']:
+            ChannelStoredM3 = self.storage_channel(option)
+            HillslopeStoredM3 = self.storage_hillslope()
+            self.var.TotalWaterStorageMM = (ChannelStoredM3 + HillslopeStoredM3) * self.var.M3toMM
 
         if (not(option['InitLisflood'])) and option['repMBTs']:
 
@@ -143,22 +166,12 @@ def dynamic(self):
             # spatially variable
 
             # This is stored:
-            ChannelStoredM3 = self.var.ChanM3.copy()
-            if option['simulateLakes']:
-                ChannelStoredM3 += self.var.LakeStorageM3Balance
-            if option['simulateReservoirs']:
-                ChannelStoredM3 += self.var.ReservoirStorageM3
-            if option['simulatePolders']:
-                ChannelStoredM3 += self.var.PolderStorageM3
-
+            ChannelStoredM3 = self.storage_channel(option)
             # ChannelStoredM3=ChanM3 + ReservoirStorageM3 + LakeStorageM3 + PolderStorageM3;
             # Water stored in channel network [m3] (including reservoirs,
             # lakes, polders)
 
-            ax_veg=self.var.SoilFraction.dims.index("vegetation")            
-            Hill1 = self.var.LZ + np.sum(self.var.SoilFraction * (self.var.CumInterception + self.var.W1 + self.var.W2 + self.var.UZ),ax_veg)                                                 
-            HillslopeStoredM3 = (self.var.WaterDepth + self.var.SnowCover + Hill1 + self.var.DirectRunoffFraction * self.var.CumInterSealed) * self.var.MMtoM * self.var.PixelArea
-
+            HillslopeStoredM3 = self.storage_hillslope()
             # Water stored at hillslope elements [m3]
             # Note that W1, W2 and TotalGroundWater are defined for the pixel's permeable fraction
             # only, which is why we need to multiply with PermeableFraction to get the volumes right (no soil moisture
@@ -184,8 +197,7 @@ def dynamic(self):
                  if not option["cropsEPIC"]: # If EPIC is active, the rice fraction initialisation is handled by EPIC (setSoilFractions in EPIC_main.py)
                     self.var.SoilFraction.values[self.var.vegetation.index('Rainfed_prescribed')] += self.var.RiceFraction
 
-                 Hill1 = self.var.LZ + (self.var.SoilFraction * (self.var.CumInterception + self.var.W1 + self.var.W2 + self.var.UZ)).sum("vegetation").values
-                 HillslopeStoredM3 = (self.var.WaterDepth + self.var.SnowCover + Hill1 + self.var.DirectRunoffFraction * self.var.CumInterSealed) * self.var.MMtoM * self.var.PixelArea                
+                 HillslopeStoredM3 = self.storage_hillslope()
                  WaterStored_nextstep =np.take(np.bincount(self.var.Catchments,weights=ChannelStoredM3),self.var.Catchments)
                  WaterStored_nextstep += np.take(np.bincount(self.var.Catchments,weights=HillslopeStoredM3),self.var.Catchments)
 
@@ -273,4 +285,4 @@ def dynamic(self):
             ones = maskinfo.in_zero() + 1.0
             numpixels = np.take(np.bincount(self.var.Catchments, weights=ones),self.var.Catchments) 
             self.var.MBErrorStorage = self.var.MBError/(self.var.WaterInit)  
-            self.var.AverageFractions = SumFractions/numpixels
+            self.var.AverageFractions = SumFractions/numpixels

src/lisfloodSettings_reference.xml

@@ -109,6 +109,7 @@ You can use builtin path variables in this template and reference to other paths
     <setoption choice="1" name="repPercolationMaps"/>
     <setoption choice="1" name="repRWS"/>
     <setoption choice="1" name="repPFMaps"/>
+    <setoption choice="1" name="repTotalWaterStorageMaps"/>
     <setoption choice="1" name="repPFForestMaps"/>
 
     <setoption choice="1" name="repTotalAbs"/>
@@ -3385,14 +3386,6 @@ soil layer 2 (lower layer) [V/V] [mm3/mm3]
 </comment>
 </textvar>
 
-<textvar name="ThetaMaps" value="$(PathOut)/th">
-<comment>
-Reported volumetric soil moisture content for  combined
-soil layer 1 and 2 [V/V] [mm3/mm3]
-</comment> 
-</textvar>
-
-
 <textvar name="UZMaps" value="$(PathOut)/uz">
 <comment>
 Reported storage in upper response box [mm]
@@ -3429,6 +3422,12 @@ Reported storage in lower response box [mm]
 </comment>
 </textvar>
 
+<textvar name="TotalWaterStorageMaps" value="$(PathOut)/tws">
+<comment>
+Reported total water storage [mm]
+</comment>
+</textvar>
+
 <textvar name="DSLRMaps" value="$(PathOut)/dslr">
 <comment>
 Reported days since last rain
@@ -3483,7 +3482,6 @@ Reported interception storage, other fraction [mm]
 </comment>
 </textvar>
 
-
 <textvar name="DSLRForestMaps" value="$(PathOut)/dslrf">
 <comment>
 Reported days since last rain
@@ -3496,6 +3494,12 @@ Reported interception storage, forest fraction [mm]
 </comment>
 </textvar>
 
+<textvar name="CumInterceptionIrrigationMaps" value="$(PathOut)/cumi">
+<comment>
+Reported interception storage, irrigation fraction [mm]
+</comment>
+</textvar>
+
 <textvar name="Theta1ForestMaps" value="$(PathOut)/thfa">
 <comment>
 Reported volumetric soil moisture content for
@@ -3517,11 +3521,17 @@ soil layer 2 (lower layer) [V/V] [mm3/mm3]
 </comment>
 </textvar>
 
-<textvar name="UZForestMaps" value="$(PathOut)/uzF">
+<textvar name="UZForestMaps" value="$(PathOut)/uzf">
 <comment>
 Reported storage in upper response box [mm]
 </comment>
 </textvar>
+
+<textvar name="CumIntSealedMaps" value="$(PathOut)/cums">
+<comment>
+Reported interception (depression) storage, direct runoff (impervious) fraction [mm]
+</comment>
+</textvar>
 </group>
 
 <group>
@@ -3790,6 +3800,12 @@ Reported upper zone outflow [mm]
 </comment>
 </textvar>
 
+<textvar name="UZOutflowIrrigationMaps" value="$(PathOut)/quzi">
+<comment>
+Reported upper zone outflow [mm]
+</comment>
+</textvar>
+
 <textvar name="GwPercUZLZForestMaps" value="$(PathOut)/uz2lzF">
 <comment>
 Reported percolation from upper to lower zone [mm]
@@ -3829,6 +3845,12 @@ soil layer 2 [V/V] [mm3/mm3]
 </comment>
 </textvar>
 
+<textvar name="UZIrrigationMaps" value="$(PathOut)/uzi">
+<comment>
+Reported storage in upper response box [mm]
+</comment>
+</textvar>
+
 </group>
 
 

tests/data/LF_ETRS89_UseCase/settings/base.xml

@@ -40,6 +40,7 @@
         <setoption choice="0" name="repsimulateReservoirs"/>
         <setoption choice="0" name="repSnowMaps"/>
         <setoption choice="0" name="repPFMaps"/>
+        <setoption choice="0" name="repTotalWaterStorageMaps"/>
 
         <setoption choice="0" name="repLZMaps"/>
 
@@ -3183,14 +3184,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT
                     </comment>
                 </textvar>
 
-                <textvar name="ThetaMaps" value="$(PathOut)/th">
-                    <comment>
-                        Reported volumetric soil moisture content for combined
-                        soil layer 1 and 2 [V/V]
-                    </comment>
-                </textvar>
-
-
                 <textvar name="UZMaps" value="$(PathOut)/uz">
                     <comment>
                         Reported storage in upper response box [mm]
@@ -3226,6 +3219,12 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT
                         Reported storage in lower response box [mm]
                     </comment>
                 </textvar>
+
+                <textvar name="TotalWaterStorageMaps" value="$(PathOut)/tws">
+                    <comment>
+                        Reported total water storage [mm]
+                    </comment>
+                </textvar>
 
                 <textvar name="DSLRMaps" value="$(PathOut)/dslr">
                     <comment>
@@ -3315,7 +3314,7 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT
                     </comment>
                 </textvar>
 
-                <textvar name="UZForestMaps" value="$(PathOut)/uzF">
+                <textvar name="UZForestMaps" value="$(PathOut)/uzf">
                     <comment>
                         Reported storage in upper response box [mm]
                     </comment>
@@ -3540,6 +3539,12 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT
                     </comment>
                 </textvar>
 
+                <textvar name="UZOutflowIrrigationMaps" value="$(PathOut)/quzi">
+                    <comment>
+                        Reported upper zone outflow [mm]
+                    </comment>
+                </textvar>
+
                 <textvar name="GwPercUZLZForestMaps" value="$(PathOut)/uz2lzF">
                     <comment>
                         Reported percolation from upper to lower zone [mm]