Feature/tws merge

.gitignore

@@ -28,8 +28,9 @@ lisflood_model.egg-info
 .vscode/
 *.ipynb
 .ipynb_checkpoints/
+tws_integration_master_diff.md
 implementation_guide_scale_offset_packing.md
 lisflood_optimization_report.md
 lisflood_optimization_report_v2.md
 # version file auto-generated by setuptools_scm (do not track)
-src/lisflood/_version.py
+src/lisflood/_version.py

docs/4_Static-Maps_reservoirs-lakes/index.md

@@ -4,6 +4,8 @@ Lakes and reservoirs can be defined as a significant volume of water, which occu
 
 The modelling of lakes and reservoirs requires three maps and a set of txt files.
 
+For the generation of total water storage output (option `repTWSMaps`), two additional maps are needed to distribute the storage mass over the extent of the lakes/reservoirs (see [Lake and reservoir extent maps](#lake-and-reservoir-extent-maps) below).
+
 
 ## Lake mask map
 
@@ -20,7 +22,7 @@ The lake mask map represents the area covered by lakes and reservoirs, it is use
 | :---| :--- | :--- | :--- |
 |Global Lakes and Wetlands Database (GLWD): <br>Large Lake Polygons (Level 1) |[GLWD leve1](https://www.worldwildlife.org/publications/global-lakes-and-wetlands-database-large-lake-polygons-level-1)|2004|Global, 1:1 to 1:3 million resolution|
 |Global Lakes and Wetlands Database (GLWD): <br>Small Lake Polygons (Level 2) |[GLWD leve1](https://www.worldwildlife.org/publications/global-lakes-and-wetlands-database-large-lake-polygons-level-2)|2004|Global, 1:1 to 1:3 million resolution|
-|Fraction of inland water| It can be prepared by using<br> the methodology explained [here](../4_Static-Maps_land-use#land-use)|NA|Global|
+|Fraction of inland water| It can be prepared by using<br> the methodology explained [here](../4_Static-Maps_land-use#land-use/index.md)|NA|Global|
 
 ### Methodology
 
@@ -146,4 +148,32 @@ Finally, lake average inflow can be retrieved from observed time series (where a
 Lake maps and tables of the European 1arcmin domain and global 3arcmin domain are based on information from [HydroLAKES](https://www.hydrosheds.org/products/hydrolakes): waterbodies classified in HydroLakes as natural lake were considered for inclusion into the lakes dataset.
 Lakes included in the European 1arcmin domain had a minimum volume of 10 hm3, a minimum lake surface area of 5 km2, a minimum upstream catchment area of 50 km2.
 Lakes included in the global 3arcmin domain had a minimum volume of 100 hm3, a minimum lake surface area of 50 km2, a minimum upstream catchment area of 250 km2.
-Lake outlet width was generally measured with GIS tools; lake average inflow was computed using OS LISFLOOD CEMS EFAS and CEMS GloFAS discharge reanalysis (GloFASv4 reanalysis for GloFASv5 tables; EFASv5 naturalized flow for EFASv6 tables).
+Lake outlet width was generally measured with GIS tools; lake average inflow was computed using OS LISFLOOD CEMS EFAS and CEMS GloFAS discharge reanalysis (GloFASv4 reanalysis for GloFASv5 tables; EFASv5 naturalized flow for EFASv6 tables).
+
+
+## Lake and reservoir extent maps
+
+The lake and reservoir modules output lake/reservoir levels at the outlet locations (output files `lakeh.nc` for locations in `lakes.nc` and `res.nc`). For the generation of total water storage maps (option `repTWSMaps`), these levels have to be converted to mass changes and spatially distributed over the extent of the lakes/reservoirs. For this, two additional maps are needed that contain the lake/reservoir ID for each pixel belonging to the area covered by the respective lake/reservoir.
+
+### General map information
+
+| Map/table name | File name; type | Units; range | Description |
+| :--- | :--- | :--- | :--- |
+| Lake extent | lake_extent.nc; <br>Type: Float32 | Units: -; <br>Range: integer ID number to identify each lake | Lake ID (ID of outflow location) for each pixel belonging to the surface extent of the lake |
+| Reservoir extent | res_extent.nc; <br>Type: Float32 | Units: -; <br>Range: integer ID number to identify each reservoir | Reservoir ID (ID of outflow location) for each pixel belonging to the surface extent of the reservoir |
+
+### Methodology
+
+The creation of the lake extent map requires:
+
+1. Lake outflow location with unique lake identifier (e.g., `lakes.nc`)
+2. External local or global dataset containing the vectorized polygons of the lake outlines (e.g., in shapefile format), such as [HydroLAKES](https://www.hydrosheds.org/products/hydrolakes), as well as the surface area of the lake
+3. Table relating the unique lake identifier of OS LISFLOOD to the lake identifier in the external data set
+
+The creation of a raster file fulfilling the criteria for a consistent lake extent map (e.g., lakes with an area smaller than the grid cell size must not be neglected; lakes assigned to only one grid cell must not be overwritten by other lakes) is not straight forward, thus a standard rasterize command is not appropriate. A dedicated script should loop over each lake defined in the outlet location file and perform the following actions:
+
+- If the total lake area is too small in comparison to the grid cell area (e.g., smaller than 5% of the cell area), the outlet location pixel is set as the lake extent.
+- Otherwise, the corresponding polygon from the shapefile is intersected with the grid cells and each grid cell containing a fraction of the lake extent is checked. If the lake covers a sufficiently large fraction (e.g., more than 7%) of the grid cell, it is assigned the lake ID, else it is neglected.
+- If no intersection fraction is large enough (e.g., if the lake is very thin but long), the outlet location pixel is set as the lake extent.
+
+The creation of the reservoir extent map can be done analogously to the lake extent map using the same methodology.

docs/5_annex_output-files/index.md

@@ -148,7 +148,9 @@ To speed up the pre-run and to prevent that results are taken from the pre-run,
 | storage in lower groundwater zone       | repLZMaps              | $mm$                       | LZMaps                           | lz      |
 | number of days since last rain          | repDSLRMaps            | $days$                     | DSLRMaps <br> DSLRForestMaps                       | dslr <br> dslF  |
 | frost index                             | repFrostIndexMaps      | $\frac{°C}{days}$          | FrostIndexMaps                                     | frost           |
-| Total Water Storage                             | repTotalWaterStorageMaps      | $mm$          | TotalWaterStorageMaps                                     | tws           |
+| Total Water Storage                             | repTWSMaps      | $m$          | TWSMaps                                     | TotalWS           |
+| Individual Water Storage Compartments           | repStorageMaps  | $m$          | LakeSMaps <br> RiverSMaps <br> SoilSMaps <br> GWSMaps <br> SnowSMaps <br> CumSMaps | Lakestor <br> Riverstor <br> Soilstor <br> GWstor <br> Snowstor <br> Cumstor          |
+| Flow Momentum                           | repFlowMomMaps  | $\frac{kgm}{s}$          | FlowMomMaps                                 | FlowMomentum  |
 | **RATE VARIABLES**                 |                        |                            |                                                    |                 |
 | rain (excluding snow)                   | repRainMaps            | $\frac{mm}{timestep}$      | RainMaps                                           | rain            |
 | snow                                    | repSnowMaps            | $\frac{mm}{timestep}$      | SnowMaps                                           | snow            |
@@ -181,8 +183,9 @@ The users should be aware that some state maps are generated only if the relevan
 
 **Note**
 
-Some cumulative storages and volumes are computed internally by LISFLOOD. Some relevant example is described below:
-- Total Water Storage is the total water volume stored in channels, lakes, reservoirs, snow cover, sealed surfaces depressions, surface runoff, canopy interception, upper and lower groundwater zones. 
+Some cumulative storages and volumes are computed internally by LISFLOOD. Some relevant examples are described below:
+- Total Water Storage is the total water mass (expressed in equivalent water height) stored in channels, lakes, reservoirs, all soil layers, snow cover, sealed surfaces depressions, surface runoff, canopy interception, uppper and lower groundwater zones. Using the option "repStorageMaps" outputs the individual storage compartments (Rivers, Lakes/Reservoirs, Soil, Groundwater, Snow, Interception storage) which add up to TWS ("repTWSMaps"). For the computation of mass stored in lakes and reservoirs (i.e. "simulateLakes" and/or "simulateReservoirs" is set to 1), the initial LISFLOOD output of lake levels centered at the outlet location are distributed over the actual extents of the lakes/reservoirs, which requires two additional input maps when using the options "repTWSMaps" and "repStorageMaps" (see [Reservoirs and lakes](../4_Static-Maps_reservoirs-lakes/)).
+- Flow Momentum is the product of flow velocity and river mass, which is needed for special applications, e.g. the computation of Earth rotation excitation functions. 
 - Surface runoff is the sum of direct runoff (from sealed and water fractions) and runoff generated by the pervious land cover fractions (forest, irrigation, other).
 - Total runoff is the sum of surface runoff and sub-surface runoff. Sub-surface runoff is the outflow from upper and lower groundwater zones.
 

src/lisflood/Lisflood_dynamic.py

@@ -248,6 +248,10 @@ def splitlanduse(array1, array2=None, array3=None):
         # Calculate water level
         self.waterlevel_module.dynamic()
 
+        # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        # Calculate water storage
+        self.waterstorage_module.dynamic()        
+
         # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
         # ************************************************************

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
@@ -208,6 +210,7 @@ def __init__(self):
         self.reservoir_module.initial()
         self.lakes_module.initial()
         self.polder_module.initial()
+        self.waterstorage_module.initial()
 
         self.transmission_module.initial()