From b38854c69b514ef4dd07388f9fa13bd3cebaf037 Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 13 Mar 2025 11:30:35 +0100 Subject: [PATCH 1/8] Update lisfloodSettings_reference.xml for TWS output --- src/lisfloodSettings_reference.xml | 62 +++++++++++++++++++++++++++++- 1 file changed, 61 insertions(+), 1 deletion(-) 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 + + + + @@ -3276,6 +3280,62 @@ Multiplier [] applied to ChanSdXdY Reported transpiration maps, weighted sum over the fractions of each pixel [mm] + + + + Reported total water storage storage [m] + added by GFZ + + + + + + Reported lake and reservoir storage [m] + added by GFZ + + + + + + Reported river storage [m] + added by GFZ + + + + + + Reported soil storage [m] + added by GFZ + + + + + + Reported groundwater storage [m] + added by GFZ + + + + + + Reported snow storage [m] + added by GFZ + + + + + + Reported cumulative interception and depression storage storage [m] + added by GFZ + + + + + + Reported riverflow momentum for HAM motion term [kgm/s] + added by GFZ + + @@ -6166,4 +6226,4 @@ Multiplier [] applied to ChanSdXdY ReportSteps=$(ReportSteps); - \ No newline at end of file + From 920e17925055c32e27338527ddecdb556a2b7b0d Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 13 Mar 2025 11:34:04 +0100 Subject: [PATCH 2/8] Update Lisflood_initial.py for TWS output --- src/lisflood/Lisflood_initial.py | 3 +++ 1 file changed, 3 insertions(+) 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() From e8e031c076a30c385008d477d119277acae41a10 Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 13 Mar 2025 11:36:28 +0100 Subject: [PATCH 3/8] Update Lisflood_dynamic.py for TWS output --- src/lisflood/Lisflood_dynamic.py | 4 ++++ 1 file changed, 4 insertions(+) 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() + # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ************************************************************ From df9dd960800c3b66aeff62b98d6ae3519b9acead Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 13 Mar 2025 11:44:25 +0100 Subject: [PATCH 4/8] Update default_options.py for TWS output --- src/lisflood/global_modules/default_options.py | 9 +++++++++ 1 file changed, 9 insertions(+) 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)}) From f085ec9ff0c439ad80f999f0b1877439a16aca1e Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 13 Mar 2025 11:47:37 +0100 Subject: [PATCH 5/8] Update routing.py for computing river flow momentum --- src/lisflood/hydrological_modules/routing.py | 5 +++++ 1 file changed, 5 insertions(+) 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 Date: Thu, 13 Mar 2025 11:54:08 +0100 Subject: [PATCH 6/8] Update evapowater.py to deal with modified LakeMask TWS computation needs a LakeMask with values other than 0 and 1. The additional line sets LakeMask back to 0/1, no impact if original (default) LakeMask is used. --- src/lisflood/hydrological_modules/evapowater.py | 1 + 1 file changed, 1 insertion(+) 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) From 1e524b75d2beba7f9ad31947ba324a1ac6119d17 Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 13 Mar 2025 11:56:47 +0100 Subject: [PATCH 7/8] Add waterstorage.py for TWS computation --- .../hydrological_modules/waterstorage.py | 284 ++++++++++++++++++ 1 file changed, 284 insertions(+) create mode 100644 src/lisflood/hydrological_modules/waterstorage.py diff --git a/src/lisflood/hydrological_modules/waterstorage.py b/src/lisflood/hydrological_modules/waterstorage.py new file mode 100644 index 00000000..e73e291d --- /dev/null +++ b/src/lisflood/hydrological_modules/waterstorage.py @@ -0,0 +1,284 @@ +""" + +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 + # gaps within lakeIDs or reservoirIDs have to be smaller than + 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_soilM = ((self.var.Theta1a[0] * self.var.SoilDepth1a[0] + self.var.Theta1b[0] * self.var.SoilDepth1b[0] + self.var.Theta2[0] * self.var.SoilDepth2[0]) * self.var.OtherFraction + + (self.var.Theta1a[1] * self.var.SoilDepth1a[1] + self.var.Theta1b[1] * self.var.SoilDepth1b[1] + self.var.Theta2[1] * self.var.SoilDepth2[1]) * self.var.ForestFraction + + (self.var.Theta1a[2] * self.var.SoilDepth1a[2] + self.var.Theta1b[2] * self.var.SoilDepth1b[2] + self.var.Theta2[2] * self.var.SoilDepth2[2]) * self.var.IrrigationFraction + ) / 1000 + + # 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 + From 6c12856e5c21eb8a528e2b6d743861dc229666b1 Mon Sep 17 00:00:00 2001 From: JensenLaura Date: Thu, 27 Mar 2025 09:57:46 +0100 Subject: [PATCH 8/8] Update waterstorage.py (soil storage computation) Soil storage is computed for each of the three soil layers separately to enable an output of individual soil layers, not only the sum of the three. --- .../hydrological_modules/waterstorage.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/lisflood/hydrological_modules/waterstorage.py b/src/lisflood/hydrological_modules/waterstorage.py index e73e291d..ac6a5c26 100644 --- a/src/lisflood/hydrological_modules/waterstorage.py +++ b/src/lisflood/hydrological_modules/waterstorage.py @@ -251,10 +251,22 @@ def dynamic(self): tws_oflowM[reservoir_mask] = np.nansum(tws_oflowM3[reservoir_mask]) / grid_area_reservoir # soil water storage [mm] -> [m] - tws_soilM = ((self.var.Theta1a[0] * self.var.SoilDepth1a[0] + self.var.Theta1b[0] * self.var.SoilDepth1b[0] + self.var.Theta2[0] * self.var.SoilDepth2[0]) * self.var.OtherFraction + - (self.var.Theta1a[1] * self.var.SoilDepth1a[1] + self.var.Theta1b[1] * self.var.SoilDepth1b[1] + self.var.Theta2[1] * self.var.SoilDepth2[1]) * self.var.ForestFraction + - (self.var.Theta1a[2] * self.var.SoilDepth1a[2] + self.var.Theta1b[2] * self.var.SoilDepth1b[2] + self.var.Theta2[2] * self.var.SoilDepth2[2]) * self.var.IrrigationFraction + 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 +