Repository: Deltares/Wflow.jl (language: python)

Issue:
Include water demand and allocation
See also https://github.com/Deltares/hydromt_wflow/issues/87
improvement reduction function Feddes
This is an issue from [wflow Python](https://github.com/openstreams/wflow/issues/24):

Apply Feddes wet reduction only to cropland (and not to natural vegetation & forest)


Discussion / hints:

Fix patch:
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index c37e167a6..252be48b4 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -70,8 +70,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   built-in operators is still allowed. This change in naming convention is now in effect and
   all invalid uses of non-ASCII characters have been replaced by ASCII equivalents.
 - Docs: 1) improved description of different model configurations in model-setup.md, also in
-  relation to hydromt_wflow in docs, 2) citing info related to wflow\_sbm publication in
+  relation to hydromt\_wflow in docs, 2) citing info related to wflow\_sbm publication in
   Geosci. Model Dev. (from in review to published).
+- Root water uptake reduction (Feddes): 1) extend critical pressure head parameters with
+  `h3_low` and `h3_high`, 2) all critical pressure head parameters can be set (values for
+  `h1`, `h2` and `h3` were fixed) and 3) the root water uptake reduction coefficient
+  ``\alpha`` can be set at 0 (default is 1) at critical pressure head `h1` (through the
+  model parameter `alpha_h1`).
+- For the actual transpiration computation in `SBM`, the potential transpiration is
+  partitioned over the soil layers with roots according to the model parameter
+  `rootfraction` (fraction of the total root length in a soil layer). Previously,
+  for each soil layer (from top to bottom) the actual transpiration was computed, and the
+  remaining potential transpiration was used in the next soil layer.
 
 ### Added
 - Total water storage as an export variable for `SBM` concept. This is the total water stored
@@ -84,6 +94,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   the following sections: [The SBM soil water accounting scheme](@ref) and [Subsurface flow
   routing](@ref) for a short description.
 - Local inertial routing to `sbm_gwf` model type.
+- Water demand and allocation computations for model types `sbm` and `sbm_gwf`.
 
 ## v0.7.3 - 2024-01-12
 
diff --git a/docs/src/images/paddy_profile.png b/docs/src/images/paddy_profile.png
new file mode 100644
index 000000000..39a134503
Binary files /dev/null and b/docs/src/images/paddy_profile.png differ
diff --git a/docs/src/index.md b/docs/src/index.md
index 3a66f391d..eba809525 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -6,13 +6,14 @@ CurrentModule = Wflow
 
 Wflow is Deltares' solution for modelling hydrological processes, allowing users to account
 for precipitation, interception, snow accumulation and melt, evapotranspiration, soil water,
-surface water and groundwater recharge in a fully distributed environment. Successfully
-applied worldwide for analyzing flood hazards, drought, climate change impacts and land use
-changes, wflow is growing to be a leader in hydrology solutions. Wflow is conceived as a
-framework, within which multiple distributed model concepts are available, which maximizes
-the use of open earth observation data, making it the hydrological model of choice for data
-scarce environments. Based on gridded topography, soil, land use and climate data, wflow
-calculates all hydrological fluxes at any given grid cell in the model at a given time step.
+surface water, groundwater recharge, and water demand and allocation in a fully distributed
+environment. Successfully applied worldwide for analyzing flood hazards, drought, climate
+change impacts and land use changes, wflow is growing to be a leader in hydrology solutions.
+Wflow is conceived as a framework, within which multiple distributed model concepts are
+available, which maximizes the use of open earth observation data, making it the
+hydrological model of choice for data scarce environments. Based on gridded topography,
+soil, land use and climate data, wflow calculates all hydrological fluxes at any given grid
+cell in the model at a given time step.
 
 Wflow was born out of the creation of Deltares in 2008, when a strategic review identified
 the need for a distributed hydrological model to allow the simulation of flows at the
diff --git a/docs/src/model_docs/lateral/kinwave.md b/docs/src/model_docs/lateral/kinwave.md
index 3e69a83f6..d0e744e73 100644
--- a/docs/src/model_docs/lateral/kinwave.md
+++ b/docs/src/model_docs/lateral/kinwave.md
@@ -54,6 +54,13 @@ External water (supply/abstraction) `inflow` [m``^3`` s``^{-1}``]  can be added
 kinematic wave for surface water routing, as a cyclic parameter or as part of forcing (see
 also [Input section](@ref)).
 
+## Abstractions
+Abstractions from the river through the variable `abstraction` [m``^3`` s``{-1}``] are
+possible when water demand and allocation is computed. The variable `abstraction` is set
+from the water demand and allocation module each time step. The `abstraction` is divided by
+the length of the runoff pathway and subtracted from the lateral inflow of the kinematic
+wave routing scheme for river flow.
+
 ## Subsurface flow routing
 In the SBM model the kinematic wave approach is used to route subsurface flow laterally.
 Different vertical hydraulic conductivity depth profiles are possible as part of the
diff --git a/docs/src/model_docs/lateral/local-inertial.md b/docs/src/model_docs/lateral/local-inertial.md
index 884e20534..ed320c754 100644
--- a/docs/src/model_docs/lateral/local-inertial.md
+++ b/docs/src/model_docs/lateral/local-inertial.md
@@ -146,6 +146,12 @@ External water (supply/abstraction) `inflow` [m``^3`` s``^{-1}``]  can be added
 inertial model for river flow (1D) and river and overland flow combined (1D-2D), as a cyclic
 parameter or as part of forcing (see also [Input section](@ref)).
 
+## Abstractions
+Abstractions from the river through the variable `abstraction` [m``^3`` s``{-1}``] are
+possible when water demand and allocation is computed. The variable `abstraction` is set
+from the water demand and allocation module each time step. Abstractions are subtracted as
+part of the continuity equation of the local inertial model.
+
 ## Multi-Threading
 The local inertial model for river flow (1D) and river and overland flow combined (1D-2D)
 can be executed in parallel using multiple threads.
diff --git a/docs/src/model_docs/model_configurations.md b/docs/src/model_docs/model_configurations.md
index f2872153c..814ddbfd2 100644
--- a/docs/src/model_docs/model_configurations.md
+++ b/docs/src/model_docs/model_configurations.md
@@ -13,6 +13,14 @@ for the vertical concept SBM of wflow\_sbm:
 - The unsaturated zone can be split-up in different layers
 - The addition of evapotranspiration losses
 - The addition of a capillary rise
+- The addition of water demand and allocation
+
+The water demand and allocation computations are supported by the wflow\_sbm model
+configurations:
+- [SBM + Kinematic wave](@ref config_sbm)
+- [SBM + Groundwater flow](@ref config_sbm_gwf)
+- [SBM + Local inertial river](@ref config_sbm_gwf_lie_river)
+- [SBM + Local inertial river (1D) and land (2D)](@ref config_sbm_gwf_lie_river_land)
 
 The vertical SBM concept is explained in more detail in the following section [SBM vertical
 concept](@ref vert_sbm).
diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index 05d89cea3..6456500e8 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -22,7 +22,9 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 | `q_av`          | average discharge | m``^3`` s``^{-1}``| - |
 | `qlat`          | lateral inflow per unit length  | m``^2`` s``^{-1}``| - |
 | `inwater`       | lateral inflow | m``^3`` s``^{-1}``| - |
-| `inflow`        | external inflow (abstraction/supply/demand) | m``^3`` s``^{-1}``| 0.0 |
+| **`inflow`**    | external inflow (abstraction/supply/demand) | m``^3`` s``^{-1}``| 0.0 |
+| `inflow_wb`     | inflow waterbody (lake or reservoir model) from land part | m``^3`` s``^{-1}``| 0.0 |
+| `abstraction`   | abstraction (computed as part of water demand and allocation) | m``^3`` s``^{-1}``| 0.0 |
 | `volume`        | kinematic wave volume |m``^3``| - |
 | `h`             | water level | m | - |
 | `h_av`          | average water level | m | - |
@@ -38,6 +40,7 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 | `lake_index`   |  map cell to 0 (no lake) or i (pick lake i in lake field) | - | - |
 | `reservoir`    | an array of reservoir models `SimpleReservoir` | - | - |
 | `lake`         | an array of lake models `Lake` | - | - |
+| `allocation`| water allocation of type `AllocationRiver` | - | - |
 | `kinwave_it`   | boolean for kinematic wave iterations | - | false |
 
 The Table below shows the parameters (fields) of struct `SurfaceFlowLand` used for overland
@@ -210,6 +213,7 @@ the `layered_exponential` profile `kv` is used and `z_exp` is required as part o
 | `ssfin` | inflow from upstream cells | m``^3`` d``{-1}``  | - |
 | `ssfmax` | maximum subsurface flow | m``^2`` d``{-1}``  | - |
 | `to_river` | part of subsurface flow that flows to the river | m``^3`` d``{-1}``  | - |
+| `volume` | subsurface water volume | m``^3`` | - |
 
 ## Local inertial
 
@@ -270,7 +274,8 @@ model parameter `mannings_n`, and is listed in the Table below between parenthes
 | `volume`    | river volume | m``^3`` | - |
 | `error`    | error volume | m``^3`` | - |
 | `inwater`    | lateral inflow | m``^3`` s``^{-1}`` | - |
-| `inflow`        | external inflow (abstraction/supply/demand) | m``^3`` s``^{-1}``| 0.0 |
+| **`inflow`**  | external inflow (abstraction/supply/demand) | m``^3`` s``^{-1}``| 0.0 |
+| `abstraction` | abstraction (computed as part of water demand and allocation) | m``^3`` s``^{-1}``| 0.0 |
 | `inflow_wb`        | inflow waterbody (lake or reservoir model) from land part | m``^3`` s``^{-1}``| 0.0 |
 | `bankfull_volume`    | bankfull volume | m``^3`` | - |
 | **`bankfull_depth`**    | bankfull depth | m | - |
@@ -280,6 +285,7 @@ model parameter `mannings_n`, and is listed in the Table below between parenthes
 | `waterbody`   |  water body cells (reservoir or lake) | - | - |
 | `reservoir`    | an array of reservoir models `SimpleReservoir` | - | - |
 | `lake` | an array of lake models `Lake` | - | - |
+| `allocation`| optional water allocation of type `AllocationRiver` | - | - |
 | `floodplain` | optional 1D floodplain routing `FloodPlain` | - | - |
 
 ### [1D floodplain](@id local-inertial_floodplain_params)
@@ -376,6 +382,18 @@ internal model parameter `z`, and is listed in the Table below between parenthes
 | `rivercells`  |  river cells| - | - |
 | `h_av` | average water depth| m | - |
 
+## Water allocation river
+The Table below shows the parameters (fields) of struct `AllocationRiver`, used when water
+demand and allocation is computed (optional), including a description of these parameters,
+the unit, and default value if applicable.
+
+|  parameter  | description  	  | unit  | default |
+|:--------------- | ------------------| ----- | -------- |
+| `act_surfacewater_abst` |  actual surface water abstraction | mm Δt⁻¹ | - |
+| `act_surfacewater_abst_vol`| actual surface water abstraction | m``^3`` Δt⁻¹ | - |
+| `available_surfacewater`| available surface water | m``^3`` | - |
+| `nonirri_returnflow`| return flow from non-irrigation (industry, domestic and livestock) | mm Δt⁻¹ | - |
+
 ## Groundwater flow
 
 ### Confined aquifer
diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md
index 02477283c..12e52bcd6 100644
--- a/docs/src/model_docs/params_vertical.md
+++ b/docs/src/model_docs/params_vertical.md
@@ -51,7 +51,7 @@ profile `kv` is used and `z_layered` is required as input.
 | **`f`** | scaling parameter (controls exponential decline of `kv_0`) | mm``^{-1}`` | 0.001  |
 | **`z_exp`** | Depth from soil surface for which exponential decline of `kv_0` is valid | mm | -  |
 | **`z_layered`** | Depth from soil surface for which layered profile (of `layered_exponential`) is valid | mm | -  |
-| **`hb`** | air entry pressure of soil (Brooks-Corey) | cm | 10.0  |
+| **`hb`** | air entry pressure of soil (Brooks-Corey) | cm | -10.0  |
 | **`soilthickness`** | soil thickness | mm | 2000.0  |
 | **`infiltcappath`** | infiltration capacity of the compacted areas  | mm Δt``^{-1}`` | 10.0 mm day``^{-1}`` |
 | **`infiltcapsoil`** | soil infiltration capacity | mm Δt``^{-1}`` | 100.0 mm day``^{-1}``|
@@ -61,6 +61,13 @@ profile `kv` is used and `z_layered` is required as input.
 | **`waterfrac`** | fraction of open water (excluding rivers) | - | 0.0  |
 | **`pathfrac`** | fraction of compacted area | - | 0.01  |
 | **`rootingdepth`** | rooting depth  | mm | 750.0  |
+| **`rootfraction`** | fraction of the root length density in each soil layer | - | - |
+| **`h1`** | soil water pressure head h1 of the root water uptake reduction function (Feddes) | cm | 0.0 cm |
+| **`h2`** | soil water pressure head h2 of the root water uptake reduction function (Feddes) | cm | -100.0 cm |
+| **`h3_high`** | soil water pressure head h3_high of the root water uptake reduction function (Feddes) | cm | -400.0 cm |
+| **`h3_low`** | soil water pressure head h3_low of the root water uptake reduction function (Feddes) | cm | -1000.0 cm |
+| **`h4`** | soil water pressure head h4 of the root water uptake reduction function (Feddes) | cm | -15849.0 cm |
+| **`alpha_h1`** |  root water uptake reduction at soil water pressure head h1 (0.0 or 1.0) | - | 1.0 |
 | **`rootdistpar`** | controls how roots are linked to water table | - | -500.0  |
 | **`cap_hmax`** | water depth beyond which capillary flux ceases | mm | 2000.0  |
 | **`cap_n`** | coefficient controlling capillary rise | - | 2.0  |
@@ -134,7 +141,12 @@ profile `kv` is used and `z_layered` is required as input.
 | `waterlevel_land` | water level land | mm | - |
 | `waterlevel_river` | water level river | mm | - |
 | `total_storage` | total water storage (excluding floodplains, lakes and reservoirs) | mm | - |
-
+| `paddy` | optional paddy (rice) fields of type `Paddy` (water demand and irrigation) | - | - |
+| `nonpaddy` | optional non-paddy fields of type `NonPaddy` (water demand and irrigation) | - | - |
+| `domestic` | optional domestic water demand of type `NonIrrigationDemand` | - | - |
+| `livestock` | optional livestock water demand of type `NonIrrigationDemand` | - | - |
+| `industry` | optional industry water demand of type `NonIrrigationDemand` | - | - |
+| `allocation` | optional water allocation of type `AllocationLand` | - | - |
 
 ## [HBV](@id params_hbv)
 The Table below shows the parameters (fields) of struct `HBV`, including a description of
@@ -378,3 +390,82 @@ specific_leaf = "Sl"
 | `TCsand` | transport capacity of overland flow for particle class sand  | ton Δt``^{-1}`` | - |
 | `TCsagg` | transport capacity of overland flow for particle class small aggregates  | ton Δt``^{-1}`` | - |
 | `TClagg` | transport capacity of overland flow for particle class large aggregates  | ton Δt``^{-1}`` | - |
+
+## Water demand and allocation
+
+### Paddy
+The Table below shows the parameters (fields) of struct `Paddy`, including a description of
+these parameters, the unit, and default value if applicable. The parameters in bold
+represent model parameters that can be set through static and forcing input data (netCDF),
+and can be listed in the TOML configuration file under `[input.vertical.paddy]`, to map the
+internal model parameter to the external netCDF variable.
+
+|  parameter | description    | unit | default |
+|:---------------| --------------- | ---------------------- | ----- |
+| `demand_gross` | irrigation gross demand | mm Δt``^{-1}`` | - |
+| **`irrigation_efficiency`** | irrigation efficiency | - | - |
+| **`maximum_irrigation_rate`** | maximum irrigation rate | mm Δt``^{-1}`` | 25.0 mm day``^{-1}`` |
+| **`irrigation_areas`** | irrigation areas | - | - |
+| **`irrigation_trigger`** | irrigation on or off (boolean) | - | - |
+| **`h_min`** | minimum required water depth in the irrigated paddy fields | mm | 20.0 |
+| **`h_opt`** | optimal water depth in the irrigated paddy fields | mm | 50.0 |
+| **`h_max`** | water depth when paddy field starts spilling water (overflow) | mm | 80.0 |
+| `h` | actual water depth in paddy field | mm | - |
+
+### Non-paddy
+The Table below shows the parameters (fields) of struct `NonPaddy`, including a description
+of these parameters, the unit, and default value if applicable. The parameters in bold
+represent model parameters that can be set through static and forcing input data (netCDF),
+and can be listed in the TOML configuration file under `[input.vertical.nonpaddy]`, to map
+the internal model parameter to the external netCDF variable.
+
+|  parameter | description    | unit | default |
+|:---------------| --------------- | ---------------------- | ----- |
+| `demand_gross` | irrigation gross demand | mm Δt``^{-1}`` | - |
+| **`irrigation_efficiency`** | irrigation efficiency | - | - |
+| **`maximum_irrigation_rate`** | maximum irrigation rate | mm Δt``^{-1}`` | 25.0 mm day``^{-1}``|
+| **`irrigation_areas`** | irrigation areas | - | - |
+| **`irrigation_trigger`** | irrigation on or off (boolean) | - | - |
+
+### Non-irrigation (industry, domestic and livestock)
+The Table below shows the parameters (fields) of struct `NonIrrigationDemand`, including a
+description of these parameters, the unit, and default value if applicable. The parameters
+in bold represent model parameters that can be set through static and forcing input data
+(netCDF). These parameters can be listed for the sectors industry, domestic and livestock,
+in the TOML configuration file under `[input.vertical.industry]`,
+`[input.vertical.domestic]` and `[input.vertical.livestock]`, to map the internal model
+parameter to the external netCDF variable.
+
+|  parameter | description    | unit | default |
+|:---------------| --------------- | ---------------------- | ----- |
+| **`demand_gross`** | gross industry water demand | mm Δt``^{-1}`` | 0.0 |
+| **`demand_net`** | net industry water demand | mm Δt``^{-1}`` | 0.0 |
+| `returnflow_fraction` | return flow fraction | - | - |
+| `returnflow` | return flow | mm Δt``^{-1}`` | - |
+
+### Water allocation land
+The Table below shows the parameters (fields) of struct `AllocationLand`, including a
+description of these parameters, the unit, and default value if applicable. The parameters
+in bold represent model parameters that can be set through static and forcing input data
+(netCDF), and can be listed in the TOML configuration file under
+`[input.vertical.allocation]`, to map the internal model parameter to the external netCDF
+variable.
+
+|  parameter | description    | unit | default |
+|:---------------| --------------- | ---------------------- | ----- |
+| `irri_demand_gross` | irrigation gross demand | mm Δt``^{-1}`` | - |
+| `nonirri_demand_gross` | non-irrigation gross demand | mm Δt``^{-1}`` | - |
+| `total_gross_demand` | total gross demand | mm Δt``^{-1}`` | - |
+| **`frac_sw_used`** | fraction surface water used | - | 1.0 |
+| **`areas`** | allocation areas | - | 1 |
+| `surfacewater_demand` | demand from surface water | mm Δt``^{-1}`` | - |
+| `surfacewater_alloc` | allocation from surface water | mm Δt``^{-1}`` | - |
+| `act_groundwater_abst` | actual groundwater abstraction | mm Δt``^{-1}`` | - |
+| `act_groundwater_abst_vol` | actual groundwater abstraction | m``^3`` Δt``^{-1}`` | - |
+| `available_groundwater` | available groundwater | m``^3`` | - |
+| `groundwater_demand` | groundwater_demand |mm Δt``^{-1}`` | - |
+| `groundwater_alloc` | allocation from groundwater |mm Δt``^{-1}`` | - |
+| `irri_alloc` | allocated water for irrigation |mm Δt``^{-1}`` | - |
+| `nonirri_alloc` | allocated water for non-irrigation |mm Δt``^{-1}`` | - |
+| `total_alloc` | total allocated water |mm Δt``^{-1}`` | - |
+| `nonirri_returnflow` | return flow from non-irrigation |mm Δt``^{-1}`` | - |
\ No newline at end of file
diff --git a/docs/src/model_docs/vertical/sbm.md b/docs/src/model_docs/vertical/sbm.md
index ca6146f76..140728566 100644
--- a/docs/src/model_docs/vertical/sbm.md
+++ b/docs/src/model_docs/vertical/sbm.md
@@ -183,48 +183,72 @@ availability.
 
 ### Transpiration
 
-The fraction of wet roots is determined using a sigmoid fuction (see figure below). The
-parameter `rootdistpar` defines the sharpness of the transition between fully wet and fully
-dry roots. The returned wet roots fraction is multiplied by the potential evaporation (and
-limited by the available water in saturated zone) to get the transpiration from the
-saturated part of the soil. This is implemented using the following code (`i` refers to the
-index of the vector that contains all active cells within the spatial model domain):
+The maximum possible root water extraction rate for each soil layer is determined by
+partitioning the potential transpiration rate ``T_p`` based on the fraction of the total
+root length (`rootfraction` [-]) in each soil layer. A root water uptake reduction model is
+used to calculate a reduction coefficient as a function of soil water pressure, that may
+reduce the maximum possible root water extraction rate. The root water uptake reduction
+model is based on the concept proposed by Feddes et al. (1978). This concept defines a
+reduction coefficient ``\alpha`` [-] as a function of soil water pressure (``h`` [cm]). Four
+different levels of ``h`` are defined: `h1`, `h2`, `h3` and `h4`. `h1` represents anoxic
+moisture conditions, `h2` represents field capacity, `h3` represents the point of critical
+soil moisture content (onset of drought stress), and `h4` represents the wilting point. The
+value of `h3` is a function of the potential transpiration rate, between 1 and 5 mm
+d``^{-1}``. If ``T_p \le 1 \text{ mm d}^{-1}``, `h3` is set equal to `h3_low` (input model
+parameter). If ``T_p \ge 5 \text{ mm d}^{-1}``, `h3` is set equal to `h3_high` (input model
+parameter). For ``T_p`` values between 1 and 5 mm d``^{-1}``, the value of `h3` is linearly
+related to ``T_p`` (between `h3_low` and `h3_high`). Besides model parameters `h3_high` and
+`h3_low`, the critical pressure heads `h1`, `h2` and `h4` can be defined as input to the
+model.
+
+The current soil water pressure is determined following the concept defined by Brooks and
+Corey (1964):
 
-```julia
-    # transpiration from saturated store
-    wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
-    actevapsat = min(sbm.pottrans[i] * wetroots, satwaterdepth)
-    satwaterdepth = satwaterdepth - actevapsat
-    restpottrans = sbm.pottrans[i] - actevapsat
+```math
+    \frac{(\theta-\theta_r)}{(\theta_s-\theta_r)} =  \Bigg\lbrace{\left(\frac{h_b}{h}\right)^{\lambda}, h > h_b \atop 1 , h \leq h_b}
 ```
 
-![soil_wetroots](../../images/soil_wetroots.png)
+where ``h`` is the pressure head [cm], ``h_b`` is the air entry pressure head [cm], and
+``\theta``, ``\theta_s``, ``\theta_r`` and ``\lambda`` as previously defined.
 
-*Amount of wet roots and the effect of the rootdistpar parameter*
+Whenever the current soil water pressure drops below `h4`, the root water uptake is set to
+zero. The root water uptake is at ideal conditions whenever the soil water pressure is above
+`h3`, with a linear transition between `h3` and `h4`. The assumption that very wet
+conditions do not affect root water uptake too much is probably generally applicable to
+natural vegetation. For crops this assumption is not valid and in this case root water
+uptake above `h1` should be set to zero (oxygen deficit) and between `h1` and `h2` root
+water uptake is limited. This is possible by setting the input model parameter `alpha_h1` at
+0 (default is 1).
+
+![soil_rootwateruptake](../../images/soil_rootwateruptake.png)
+
+*Root water uptake reduction coefficient as a function of soil water pressure*
 
 ```@setup
 # Figure created using python: # hide
 # https://gist.github.com/JoostBuitink/21dd32e71fd1360117fcd1c532c4fd9d#file-sbm_soil_figs-py # hide
 ```
 
-The remaining potential evaporation is used to extract water from the unsaturated store. The
-maximum allowed extraction of the unsaturated zone is determined based on the fraction of
-the unsaturated zone that is above the rooting depth, see conceptual figure below. This is
-implemented using the following code:
+The maximum allowed root water extraction from each soil layer in the unsaturated zone is
+determined based on the fraction of each soil layer in the unsaturated zone that is above
+the rooting depth (`availcap`) and the unsaturated storage `usld`, see conceptual figure
+below. This is implemented using the following code (`i` refers to the index of the vector
+that contains all active cells within the spatial model domain and `k` refers to the soil
+layer (from top to bottom) in the unsaturated zone):
 
 ```julia
-    if ust # whole_ust_available = true
-        availcap = ustorelayerdepth * 0.99
+    # availcap is fraction of soil layer containing roots
+    # if `ust` is `true`, the whole unsaturated store is available for transpiration
+    if ust
+        availcap = usld[k] * 0.99
     else
-        if usl > 0.0
-            availcap = min(1.0, max(0.0, (rootingdepth - sumlayer) / usl))
-        else
-            availcap = 0.0
-        end
+        availcap =
+            min(1.0, max(0.0, (sbm.rootingdepth[i] - sbm.sumlayers[i][k]) / usl[k]))
     end
-    maxextr = availcap * ustorelayerdepth
+    maxextr = usld[k] * availcap
 ```
 
+
 ![soil_unsatevap](../../images/soil_unsatevap.png)
 
 *Conceptual overview of how maxextr depends on rooting depth and water table depth*
@@ -244,34 +268,48 @@ implemented using the following code:
     whole_ust_available = true
     ```
 
-Next, a root water uptake reduction model is used to calculate a reduction coefficient as a
-function of soil water pressure. This concept is based on the concept presented by Feddes et
-al. (1978). This concept defines a reduction coefficient `a` as a function of soil water
-pressure (`h`). Four different levels of `h` are defined: `h2`, `h3`, and `h4` are defined
-as fixed values, and `h1` can be defined as input to the model (defaults to -10 cm). `h1`
-represents the air entry pressure, `h2` represents field capacity, `h3` represents the point
-of critical soil moisture content, and `h4` represents the wilting point. The current soil
-water pressure is determined following the concept defined by Brooks and Corey (1964):
+The computation of transpiration from the saturated store depends on the water table depth,
+rooting depth, the reduction coefficient ``\alpha``, the fraction of wet roots and the
+`rootfraction` below the water table. The fraction of wet roots is determined using a
+sigmoid fuction (see figure below). The parameter `rootdistpar` defines the sharpness of the
+transition between fully wet and fully dry roots. If the water table depth is equal to or
+lower than the rooting depth, the remaining potential transpiration is used based on the
+potential transpiration and actual transpiration in the unsaturated zone. The remaining
+potential transpiration is multiplied by the wet roots fraction and the reduction
+coefficient (and limited by the available water in saturated zone) to get the transpiration
+from the saturated part of the soil. If the water table depth intersects the rooting depth,
+the potential transpiration is multiplied by the remaining `rootfraction` (below the water
+table), wet roots fraction and the reduction coefficient (and limited by the available water
+in saturated zone) to get the transpiration from the saturated part of the soil. This is
+implemented using the following code (`i` refers to the index of the vector that contains
+all active cells within the spatial model domain):
 
-```math
-    \frac{(\theta-\theta_r)}{(\theta_s-\theta_r)} =  \Bigg\lbrace{\left(\frac{h_b}{h}\right)^{\lambda}, h > h_b \atop 1 , h \leq h_b}
+```julia
+        # transpiration from saturated store
+        wetroots = scurve(sbm.zi[i], sbm.rootingdepth[i], Float(1.0), sbm.rootdistpar[i])
+        alpha = rwu_reduction_feddes(
+            Float(0.0),
+            sbm.h1[i],
+            sbm.h2[i],
+            sbm.h3[i],
+            sbm.h4[i],
+            sbm.alpha_h1[i],
+        )
+        # include remaining root fraction if rooting depth is below water table zi
+        if sbm.zi[i] >= sbm.rootingdepth[i]
+            f_roots = wetroots
+            restevap = sbm.pottrans[i] - actevapustore
+        else
+            f_roots = wetroots * (1.0 - rootfraction_unsat)
+            restevap = sbm.pottrans[i]
+        end
+        actevapsat = min(restevap * f_roots * alpha, satwaterdepth)
+        satwaterdepth = satwaterdepth - actevapsat
 ```
 
-where ``h`` is the pressure head [cm], ``h_b`` is the air entry pressure head [cm], and
-``\theta``, ``\theta_s``, ``\theta_r`` and ``\lambda`` as previously defined.
-
-Whenever the current soil water pressure drops below `h4`, the root water uptake is set to
-zero. The root water uptake is at ideal conditions whenever the soil water pressure is above
-`h3`, with a linear transition between `h3` and `h4`. Note that in the original
-transpiration reduction-curve of Feddes (1978) root water uptake above `h1` is set to zero
-(oxygen deficit) and between `h1` and `h2` root water uptake is limited. The assumption that
-very wet conditions do not affect root water uptake too much is probably generally
-applicable to natural vegetation, however for crops this assumption is not valid. This could
-be improved in the wflow code by applying the reduction to crops only.
-
-![soil_rootwateruptake](../../images/soil_rootwateruptake.png)
+![soil_wetroots](../../images/soil_wetroots.png)
 
-*Root water uptake reduction coefficient as a function of soil water pressure*
+*Amount of wet roots and the effect of the rootdistpar parameter*
 
 ```@setup
 # Figure created using python: # hide
@@ -441,15 +479,14 @@ have different infiltration capacities. Naturally, only the water that can be st
 soil can infiltrate. If not all water can infiltrate, this is added as excess water to the
 runoff routing scheme.
 
-The infiltrating
-water is split in two parts, the part that falls on compacted areas and the part that falls
-on non-compacted areas. The maximum amount of water that can infiltrate in these areas is
-calculated by taking the minimum of the maximum infiltration rate (`infiltcapsoil` [mm
-t``^{-1}``] for non-compacted areas and `infiltcappath` [mm t``^{-1}``] for compacted areas)
-and the amount of water available for infiltration `avail_forinfilt` [mm t``^{-1}``]. The
-water that can actually infiltrate `infiltsoilpath` [mm t``^{-1}``] is calculated by taking
-the minimum of the total maximum infiltration rate (compacted and non-compacted areas) and
-the remaining storage capacity.
+The infiltrating water is split in two parts, the part that falls on compacted areas and the
+part that falls on non-compacted areas. The maximum amount of water that can infiltrate in
+these areas is calculated by taking the minimum of the maximum infiltration rate
+(`infiltcapsoil` [mm t``^{-1}``] for non-compacted areas and `infiltcappath` [mm t``^{-1}``]
+for compacted areas) and the amount of water available for infiltration `avail_forinfilt`
+[mm t``^{-1}``]. The water that can actually infiltrate `infiltsoilpath` [mm t``^{-1}``] is
+calculated by taking the minimum of the total maximum infiltration rate (compacted and
+non-compacted areas) and the remaining storage capacity.
 
 Infiltration excess occurs when the infiltration capacity is smaller then the throughfall
 and stemflow rate. This amount of water (`infiltexcess` [mm t``^{-1}``]) becomes overland
@@ -587,6 +624,187 @@ Part of the water available for infiltration is diverted to the open water, base
 fractions of river and lakes of each grid cell. The amount of evaporation from open water is
 assumed to be equal to potential evaporation (if sufficient water is available).
 
+## Non-irrigation
+Non-irrigation water demand and allocation computations are supported for the sectors
+domestic, industry and livestock. These computations can be enabled by specifying the
+following in the TOML file:
+
+```toml
+[model.water_demand]
+domestic = true
+industry = true
+livestock = true
+```
+
+For these non-irrigation sectors the gross demand (``d_\mathrm{gross}`` [mm t``^{-1}``]) and
+net demand (``d_\mathrm{net}`` [mm t``^{-1}``]) are provided to the model (input through
+cyclic or forcing data). Gross demand represents the total demand and hence the total
+abstraction from surface water or groundwater when sufficient water is available. Net demand
+represents water consumption. The portion of total abstracted water that is not consumed is
+returned as surface water. The return flow fraction (``f_\mathrm{return}`` [-]) is
+calculated as follows:
+
+```math
+    f_\mathrm{return} = 1.0 - \frac{d_\mathrm{net}}{d_\mathrm{gross}},
+```
+and used to calculate the return flow rate (water abstracted from surface water or
+groundwater but not consumed). For grid cells containing a river the return flow is directly
+returned to the river routing component, otherwise the return flow is returned to the
+overland flow routing component.
+
+## Non-paddy irrigation
+Non-paddy (other crops than flooded rice) water demand and allocation computations are
+supported. These computations can be enabled by specifying the following in the TOML file:
+
+```toml
+[model.water_demand]
+nonpaddy = true
+```
+Irrigation is applied during the growing season (when input parameter `irrigation_trigger`
+[-] is `true` (or `on`)) and when water depletion exceeds the readily available water:
+
+```math
+    (U_\mathrm{field} - U_\mathrm{a}) \ge (U_\mathrm{field} - U_\mathrm{h3})
+```
+where ``U_\mathrm{field}`` \[mm\] is the unsaturated store in the root zone at field
+capacity (defined at a soil water pressure head of -100 cm), ``U_\mathrm{a}`` \[mm\] is the
+actual unsaturated store in the root zone and ``U_\mathrm{h3}`` \[mm\] is the unsaturated
+store in the root zone at the critical soil water pressure head `h3`, below this pressure
+head reduction of root water uptake starts due to drought stress. The net irrigation demand
+[mm t``^{-1}``] is the irrigation rate that brings the root zone back to field capacity,
+limited by the soil infiltration capacity [mm t``^{-1}``], assuming that farmers do not
+apply an irrigation rate higher than the soil infiltration capacity. To account for limited
+irrigation efficiency the net irrigation demand is divided by the irrigation efficiency for
+non-paddy crops (`irrigation_efficiency` [-], default is 1.0), resulting in gross irrigation
+demand [mm t``^{-1}``]. Finally, the gross irrigation demand is limited by the maximum
+irrigation rate (`maximum_irrigation_rate` [mm t``^{-1}``], default is 25 mm d``^{-1}``). If
+the maximum irrigation rate is applied, irrigation continues at subsequent time steps until
+field capacity is reached. Irrigation is added to the `SBM` variable `avail_forinfilt` [mm
+t``^{-1}``], the amount of water available for infiltration.
+
+## Paddy irrigation
+Paddy (flooded rice) water demand and allocation computations are supported. These
+computations can be enabled by specifying the following in the TOML file:
+
+```toml
+[model.water_demand]
+paddy = true
+```
+Irrigation is applied during the growing season (when input parameter `irrigation_trigger`
+[-] is `true` (or `on`)) and when the paddy water depth `h` \[mm\] reaches below the minimum
+water depth `h_min` \[mm\] (see also the figure below). The net irrigation demand [mm
+t``^{-1}``] is the irrigation rate required to reach the optimal paddy water depth `h_opt`
+\[mm\], an approach similar to Xie and Cui (2011). To account for limited irrigation
+efficiency the net irrigation demand is divided by the irrigation efficiency for paddy
+fields (`irrigation_efficiency` [-], default is 1.0), resulting in gross irrigation demand
+[mm t``^{-1}``]. Finally, the gross irrigation demand is limited by the maximum irrigation
+rate (`maximum_irrigation_rate` [mm t``^{-1}``], default is 25 mm d``^{-1}``). If the
+maximum irrigation rate is applied, irrigation continues at subsequent time steps until the
+optimal paddy water depth `h_opt` is reached. Irrigation is added to the `SBM` variable
+`avail_forinfilt` [mm t``^{-1}``], the amount of water available for infiltration. When the
+paddy water depth `h` exceeds `h_max` \[mm\] runoff occurs, and this amount is added to the
+runoff routing scheme for overland flow. The figure below shows a typical vertical soil
+profile of a puddled rice soil with a muddy layer of about 15 cm (in this case represented
+by two soil layers of 5 cm and 10 cm thickness), a plow soil layer of 5 cm with relative low
+permeability (vertical hydraulic conductivity ``k_v`` of about 5 mm d``^{-1}``), and a
+non-puddled soil below the plow soil layer. The low vertical hydraulic conductivity of the
+plow soil layer can be realized by making use of the parameter `kvfrac` [-], a
+multiplication factor applied to the vertical hydraulic conductivity at soil depth ``z``
+[mm].
+
+![paddy_profile](../../images/paddy_profile.png)
+
+*Schematic diagram of a paddy field with water balance components and soil profile*
+
+## Water withdrawal and allocation
+For the water withdrawal the total gross demand is computed (sum over the irrigation and
+non-irrigation water demand sectors), in case sufficient water is available the water
+withdrawal is equal to the total gross demand. In case of insufficient water availability,
+the water withdrawal is scaled down to the available water, and allocation is then
+proportional to the gross demand per sector (industry, domestic, livestock and irrigation).
+Water can be abstracted from the following sources:
+
+- surface water from rivers (max 80% of total available water) 
+- reservoirs and lakes (max 98% of total available water)
+- groundwater (max 75% of total available water)
+
+The model parameter `frac_sw_used` (fraction surface water used, default is 1.0) determines
+how much water is supplied by available surface water and groundwater.
+
+### Local
+First, surface water abstraction (excluding reservoir and lake locations) is computed to
+satisfy local (same grid cell) water demand. The available surface water volume is limited
+by a fixed scaling factor of 0.8 to prevent rivers from completely drying out. It is assumed
+that the water demand cannot be satisfied completely from local surface water and
+groundwater. The next step is to satisfy the remaining water demand for allocation `areas`
+[-], described in the next sub-section.
+
+### Allocation areas
+For allocation areas the water demand ``V_\mathrm{sw, demand}`` [m``^3``] and availability
+``V_\mathrm{sw, availabilty}`` [m``^3``] are summed (including reservoir and lake locations
+limited by a fixed scaling factor of 0.98), and the total surface water abstraction is then:
+
+```math
+    V_\mathrm{sw, abstraction} = \mathrm{min}(V_\mathrm{sw, demand}, V_\mathrm{sw, availabilty})
+```
+The fraction of available surface water that can be abstracted ``f_\mathrm{sw,
+abstraction}`` [-] at the allocation area level is then:
+
+```math
+    f_\mathrm{sw, abstraction} = \frac{V_\mathrm{sw, abstraction}}{V_\mathrm{sw, available}}
+```
+This fraction is applied to the remaining available surface water of each river cell
+(including lake and reservoir locations) to compute surface water abstraction at each river
+cell and to update the local surface water abstraction.
+
+The fraction of water demand that can be satisfied by available surface water
+``f_\mathrm{sw, allocation}`` [-] at the allocation area level is then:
+
+```math
+    f_\mathrm{sw, allocation} = \frac{V_\mathrm{sw, abstraction}}{V_\mathrm{sw, demand}}
+```
+This fraction is applied to the remaining surface water demand of each land cell to compute
+the allocated surface water to each land cell.
+
+Then groundwater abstraction is computed to satisfy the remaining local water demand, where
+groundwater abstraction is limited by a fixed scaling factor of 0.75 applied to the
+groundwater volume. Finally, for allocation `areas` the water demand ``V_\mathrm{gw,
+demand}`` [m``^3``] and availability ``V_\mathrm{gw, availabilty}`` [m``^3``] are summed,
+and the total groundwater abstraction is then:
+
+```math
+    V_\mathrm{gw, abstraction} = \mathrm{min}(V_\mathrm{gw, demand}, V_\mathrm{gw, availabilty})
+```
+The fraction of available groundwater that can be abstracted at allocation area level
+``f_\mathrm{gw, abstraction}`` [-] at the allocation area level is then:
+
+```math
+    f_\mathrm{gw, abstraction} = \frac{V_\mathrm{gw, abstraction}}{V_\mathrm{gw, available}}
+```
+This fraction is applied to the remaining available groundwater of each land cell to compute
+groundwater abstraction and to update the local groundwater abstraction.
+
+The fraction of water demand that can be satisfied by available groundwater ``f_\mathrm{gw,
+allocation}`` [-] at the allocation area level is then:
+
+```math
+    f_\mathrm{gw, allocation} = \frac{V_\mathrm{gw, abstraction}}{V_\mathrm{gw, demand}}
+```
+This fraction is applied to the remaining groundwater demand of each land cell to compute
+the allocated groundwater to each land cell.
+
+### Abstractions
+Groundwater abstraction is implemented by subtracting this amount from the `recharge`
+variable of the lateral subsurface flow component (kinematic wave) or the recharge `rate` of
+the groundwater flow module. Surface water `abstraction` [m``^3`` s``^{-1}``] is divided by
+the flow length `dl` [m] and subtracted from the lateral inflow of kinematic wave routing
+scheme for river flow. For the local inertial routing scheme (river and optional floodplain
+routing), the surface water `abstraction` [m``^3`` s``^{-1}``] is subtracted as part of the
+continuity equation of the local inertial model. For reservoir and lake locations surface
+water is abstracted (`act_surfacewater_abst_vol` [m``^3`` t``^{-1}``]) from the reservoir
+`volume` [m``^3``] and lake `storage` [m``^3``] respectively, with a subsequent update of
+the lake `waterlevel` [m].
+
 ## References
 + Brooks, R. H., and Corey, A. T., 1964, Hydraulic properties of porous media, Hydrology
   Papers 3, Colorado State University, Fort Collins, 27 p.
@@ -609,3 +827,5 @@ assumed to be equal to potential evaporation (if sufficient water is available).
 + Wigmosta, M. S., Lane, L. J., Tagestad, J. D., and Coleman A. M., 2009, Hydrologic and
   erosion models to assess land use and management practices affecting soil erosion, J.
   Hydrol. Eng., 14, 27-41.
++ Xie, X. and Cui, Y., 2011, Development and test of SWAT for modeling hydrological
+  processes in irrigation districts with paddy rice, J. Hydrol., 396, pp. 61-71.
diff --git a/docs/src/user_guide/additional_options.md b/docs/src/user_guide/additional_options.md
index 740943c42..6b5baa03b 100644
--- a/docs/src/user_guide/additional_options.md
+++ b/docs/src/user_guide/additional_options.md
@@ -134,6 +134,55 @@ q = "q_floodplain"
 h = "h_floodplain"
 ```
 
+## Enabling water demand and allocation
+The model types `sbm` and `sbm_gwf` support water demand and allocation computations, in
+combination with the kinematic wave and local inertial runoff routing scheme for river and
+overland flow.
+
+```toml
+# example of water demand and allocation input parameters as cyclic data
+[input]
+cyclic = ["vertical.domestic.demand_gross", "vertical.domestic.demand_net", 
+"vertical.industry.demand_gross", "vertical.industry.demand_net", 
+"vertical.livestock.demand_gross", "vertical.livestock.demand_net", 
+"vertical.paddy.irrigation_trigger", "vertical.nonpaddy.irrigation_trigger",]
+
+[model.water_demand]
+domestic = true     # optional, default is "false"
+industry = true     # optional, default is "false"
+livestock = true    # optional, default is "false"
+paddy = true        # optional, default is "false"
+nonpaddy = true     # optional, default is "false"
+
+[input.vertical.allocation]
+areas = "allocation_areas"
+frac_sw_used = "SurfaceWaterFrac"
+
+[input.vertical.domestic]
+demand_gross = "dom_gross"
+demand_net = "dom_net"
+
+[input.vertical.industry]
+demand_gross = "ind_gross"
+demand_net = "ind_net"
+
+[input.vertical.livestock]
+demand_gross = "lsk_gross"
+demand_net = "lsk_net"
+
+[input.vertical.paddy]
+irrigation_areas = "paddy_irrigation_areas"
+irrigation_trigger = "irrigation_trigger"
+
+[input.vertical.nonpaddy]
+irrigation_areas = "nonpaddy_irrigation_areas"
+irrigation_trigger = "irrigation_trigger"
+
+# required if paddy is set to "true"
+[state.vertical.paddy]
+h = "h_paddy"
+```
+
 ## [Using multithreading] (@id multi_threading)
 
 ### Using wflow in Julia
diff --git a/docs/src/user_guide/sample_data.md b/docs/src/user_guide/sample_data.md
index e35ded8b8..dff6f10f5 100644
--- a/docs/src/user_guide/sample_data.md
+++ b/docs/src/user_guide/sample_data.md
@@ -21,7 +21,7 @@ and [Command Line Interface](@ref cli).
 ```julia
 # urls to TOML and netCDF of the Moselle example model
 toml_url = "https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sbm_config.toml"
-staticmaps = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.7/staticmaps-moselle.nc"
+staticmaps = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.9/staticmaps-moselle.nc"
 forcing = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.6/forcing-moselle.nc"
 instates = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.6/instates-moselle.nc"
 
@@ -41,7 +41,7 @@ download(toml_url, toml_path)
 ```julia
 # urls to TOML and netCDF of the Moselle example model
 toml_url = "https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sbm_gwf_config.toml"
-staticmaps = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.2/staticmaps-sbm-groundwater.nc"
+staticmaps = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.3/staticmaps-sbm-groundwater.nc"
 forcing = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.1/forcing-sbm-groundwater.nc"
 
 # create a "data" directory in the current directory
diff --git a/server/Manifest.toml b/server/Manifest.toml
index d8b6f626a..aba619527 100644
--- a/server/Manifest.toml
+++ b/server/Manifest.toml
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.2"
 manifest_format = "2.0"
 project_hash = "b955683a8ff3d663c1ae626b1dc754fbae9fa2f7"
 
@@ -113,7 +113,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.0+0"
 
 [[deps.CpuId]]
 deps = ["Markdown"]
@@ -392,7 +392,7 @@ weakdeps = ["Adapt"]
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
diff --git a/src/Wflow.jl b/src/Wflow.jl
index 15784641d..f4a8bd8f5 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -114,10 +114,11 @@ struct SedimentModel end    # "sediment" type / sediment_model.jl
 Base.show(io::IO, m::Model) = print(io, "model of type ", typeof(m))
 
 include("horizontal_process.jl")
+include("flow.jl")
 include("hbv.jl")
+include("water_demand.jl")
 include("sbm.jl")
 include("flextopo.jl")
-include("flow.jl")
 include("sediment.jl")
 include("reservoir_lake.jl")
 include("hbv_model.jl")
diff --git a/src/flextopo_model.jl b/src/flextopo_model.jl
index 5575facba..039759940 100644
--- a/src/flextopo_model.jl
+++ b/src/flextopo_model.jl
@@ -635,8 +635,7 @@ function initialize_flextopo_model(config::Config)
         order = toposort,
         indices = inds,
         reverse_indices = rev_inds,
-        xl = xl,
-        yl = yl,
+        area = xl .* yl,
     )
     river = (
         graph = graph_riv,
diff --git a/src/flow.jl b/src/flow.jl
index 31fa480d7..c778b9f6c 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -1,7 +1,7 @@
 
 abstract type SurfaceFlow end
 
-@get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowRiver{T,R,L} <:
+@get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowRiver{T,R,L,W} <:
                                                                SurfaceFlow
     beta::T | "-" | 0 | "scalar"                 # constant in Manning's equation
     sl::Vector{T} | "m m-1"                      # Slope [m m⁻¹]
@@ -14,6 +14,7 @@ abstract type SurfaceFlow end
     inwater::Vector{T} | "m3 s-1"                # Lateral inflow [m³ s⁻¹]
     inflow::Vector{T} | "m3 s-1"                 # External inflow (abstraction/supply/demand) [m³ s⁻¹]
     inflow_wb::Vector{T} | "m3 s-1"              # inflow waterbody (lake or reservoir model) from land part [m³ s⁻¹]
+    abstraction::Vector{T} | "m3 s-1"            # Abstraction (computed as part of water demand and allocation) [m³ s⁻¹]
     volume::Vector{T} | "m3"                     # Kinematic wave volume [m³] (based on water level h)
     h::Vector{T} | "m"                           # Water level [m]
     h_av::Vector{T} | "m"                        # Average water level [m]
@@ -29,6 +30,7 @@ abstract type SurfaceFlow end
     lake_index::Vector{Int} | "-"                # map cell to 0 (no lake) or i (pick lake i in lake field)
     reservoir::R | "-" | 0                       # Reservoir model struct of arrays
     lake::L | "-" | 0                            # Lake model struct of arrays
+    allocation::W | "-" | 0                      # Water allocation
     kinwave_it::Bool | "-" | 0 | "none" | "none" # Boolean for iterations kinematic wave
 
     # TODO unclear why this causes a MethodError
@@ -146,6 +148,8 @@ function initialize_surfaceflow_river(
     )
     clamp!(sl, 0.00001, Inf)
 
+    do_water_demand = haskey(config.model, "water_demand")
+
     n = length(inds)
 
     sf_river = SurfaceFlowRiver(
@@ -159,6 +163,7 @@ function initialize_surfaceflow_river(
         qlat = zeros(Float, n),
         inwater = zeros(Float, n),
         inflow = zeros(Float, n),
+        abstraction = zeros(Float, n),
         inflow_wb = zeros(Float, n),
         volume = zeros(Float, n),
         h = zeros(Float, n),
@@ -176,6 +181,7 @@ function initialize_surfaceflow_river(
         reservoir = reservoir,
         lake = lake,
         kinwave_it = iterate,
+        allocation = do_water_demand ? initialize_allocation_river(n) : nothing,
     )
 
     return sf_river
@@ -288,11 +294,15 @@ function update(sf::SurfaceFlowRiver, network, doy)
                     # Inflow supply/abstraction is added to qlat (divide by flow length)
                     # If inflow < 0, abstraction is limited
                     if sf.inflow[v] < 0.0
-                        max_abstract = min(sf.qin[v] + sf.qlat[v] * sf.dl[v], -sf.inflow[v])
+                        max_abstract = min(
+                            (sf.inwater[v] + sf.qin[v] + sf.volume[v] / dt) * 0.80,
+                            -sf.inflow[v],
+                        )
                         inflow = -max_abstract / sf.dl[v]
                     else
                         inflow = sf.inflow[v] / sf.dl[v]
                     end
+                    inflow -= sf.abstraction[v] / sf.dl[v]
 
                     sf.q[v] = kinematic_wave(
                         sf.qin[v],
@@ -350,6 +360,8 @@ function update(sf::SurfaceFlowRiver, network, doy)
                     # update h
                     crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
                     sf.h[v] = crossarea / sf.width[v]
+                    sf.volume[v] = sf.dl[v] * sf.width[v] * sf.h[v]
+
                     sf.q_av[v] += sf.q[v]
                     sf.h_av[v] += sf.h[v]
                 end
@@ -358,7 +370,6 @@ function update(sf::SurfaceFlowRiver, network, doy)
     end
     sf.q_av ./= its
     sf.h_av ./= its
-    sf.volume .= sf.dl .* sf.width .* sf.h
 end
 
 function stable_timestep(sf::S) where {S<:SurfaceFlow}
@@ -372,7 +383,8 @@ function stable_timestep(sf::S) where {S<:SurfaceFlow}
             courant = zeros(n)
             for v = 1:n
                 if sf.q[v] > 0.0
-                    sf.cel[v] = 1.0 / (sf.alpha[v] * sf.beta * pow(sf.q[v], (sf.beta - 1.0)))
+                    sf.cel[v] =
+                        1.0 / (sf.alpha[v] * sf.beta * pow(sf.q[v], (sf.beta - 1.0)))
                     courant[v] = (sf.dt / sf.dl[v]) * sf.cel[v]
                 end
             end
@@ -408,6 +420,7 @@ end
     ssfin::Vector{T} | "m3 d-1"            # Inflow from upstream cells [m³ d⁻¹]
     ssfmax::Vector{T} | "m2 d-1"           # Maximum subsurface flow [m² d⁻¹]
     to_river::Vector{T} | "m3 d-1"         # Part of subsurface flow [m³ d⁻¹] that flows to the river
+    volume::Vector{T} | "m3"               # Subsurface volume [m³]
 
     function LateralSSF{T}(args...) where {T}
         equal_size_vectors(args)
@@ -415,9 +428,10 @@ end
     end
 end
 
-
 function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
-    @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network
+    @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes, area =
+        network
+
 
     ns = length(subdomain_order)
     for k = 1:ns
@@ -474,6 +488,10 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
                         ssf.ssfmax[v],
                     )
                 end
+                ssf.volume[v] =
+                    (ssf.theta_s[v] - ssf.theta_r[v]) *
+                    (ssf.soilthickness[v] - ssf.zi[v]) *
+                    area[v]
             end
         end
     end
@@ -487,7 +505,7 @@ end
     ssf::Vector{T} | "m3 d-1"           # Subsurface flow [m³ d⁻¹]
 end
 
-@get_units @exchange @grid_type @grid_location @with_kw struct ShallowWaterRiver{T,R,L,F}
+@get_units @exchange @grid_type @grid_location @with_kw struct ShallowWaterRiver{T,R,L,F,W}
     n::Int | "-" | 0 | "none" | "none"                      # number of cells
     ne::Int | "-" | 0 | "none" | "none"                     # number of edges/links
     active_n::Vector{Int} | "-"                             # active nodes
@@ -519,6 +537,7 @@ end
     error::Vector{T} | "m3"                                 # error volume
     inwater::Vector{T} | "m3 s-1"                           # lateral inflow [m³ s⁻¹]
     inflow::Vector{T} | "m3 s-1"                            # external inflow (abstraction/supply/demand) [m³ s⁻¹]
+    abstraction::Vector{T} | "m3 s-1"                       # abstraction (computed as part of water demand and allocation) [m³ s⁻¹]
     inflow_wb::Vector{T} | "m3 s-1"                         # inflow waterbody (lake or reservoir model) from land part [m³ s⁻¹]
     bankfull_volume::Vector{T} | "m3"                       # bankfull volume
     bankfull_depth::Vector{T} | "m"                         # bankfull depth
@@ -530,6 +549,7 @@ end
     reservoir::R | "-" | 0                                  # Reservoir model struct of arrays
     lake::L | "-" | 0                                       # Lake model struct of arrays
     floodplain::F | "-" | 0                                 # Floodplain (1D) schematization
+    allocation::W | "-" | 0                                 # Water allocation
 end
 
 function initialize_shallowwater_river(
@@ -668,6 +688,7 @@ function initialize_shallowwater_river(
     waterbody = !=(0).(reservoir_index .+ lake_index)
     active_index = findall(x -> x == 0, waterbody)
 
+    do_water_demand = haskey(config.model, "water_demand")
     sw_river = ShallowWaterRiver(
         n = n,
         ne = _ne,
@@ -697,6 +718,7 @@ function initialize_shallowwater_river(
         volume = fill(0.0, n),
         error = zeros(n),
         inflow = zeros(n),
+        abstraction = zeros(n),
         inflow_wb = zeros(n),
         inwater = zeros(n),
         dl = dl,
@@ -711,6 +733,7 @@ function initialize_shallowwater_river(
         reservoir = reservoir,
         lake = lake,
         floodplain = floodplain,
+        allocation = do_water_demand ? initialize_allocation_river(n) : nothing,
     )
     return sw_river, nodes_at_link
 end
@@ -885,7 +908,8 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
 
             q_src = sum_at(sw.q, links_at_node.src[i])
             q_dst = sum_at(sw.q, links_at_node.dst[i])
-            sw.volume[i] = sw.volume[i] + (q_src - q_dst + sw.inwater[i]) * dt
+            sw.volume[i] =
+                sw.volume[i] + (q_src - q_dst + sw.inwater[i] - sw.abstraction[i]) * dt
 
             if sw.volume[i] < 0.0
                 sw.error[i] = sw.error[i] + abs(sw.volume[i])
@@ -1153,8 +1177,8 @@ function stable_timestep(sw::ShallowWaterLand{T})::T where {T}
     dt_min = T(Inf)
     @batch per = thread reduction = ((min, dt_min),) for i = 1:sw.n
         @fastmath @inbounds dt =
-            sw.rivercells[i] == 0 ? sw.alpha * min(sw.xl[i], sw.yl[i]) / sqrt(sw.g * sw.h[i]) :
-            T(Inf)
+            sw.rivercells[i] == 0 ?
+            sw.alpha * min(sw.xl[i], sw.yl[i]) / sqrt(sw.g * sw.h[i]) : T(Inf)
         dt_min = min(dt, dt_min)
     end
     dt_min = isinf(dt_min) ? T(10.0) : dt_min
@@ -1323,7 +1347,7 @@ function shallowwater_update(
                         sum_at(swr.q, links_at_node.dst[inds_riv[i]]) + sw.qx[xd] -
                         sw.qx[i] + sw.qy[yd] - sw.qy[i] +
                         swr.inflow[inds_riv[i]] +
-                        sw.runoff[i]
+                        sw.runoff[i] - swr.abstraction[inds_riv[i]]
                     ) * dt
                 if sw.volume[i] < T(0.0)
                     sw.error[i] = sw.error[i] + abs(sw.volume[i])
@@ -1604,28 +1628,37 @@ function initialize_floodplain_1d(
 end
 
 """
-    set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T}
+    set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
 
-Set `inwater` of the lateral river component for a `Model` with vertical `SBM` concept.
+Set `inwater` of the lateral river component for a `Model` of type `SbmModel` or `SbmGwfModel`.
 `ssf_toriver` is the subsurface flow to the river.
 """
-function set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T}
-    @unpack lateral, vertical, network = model
+function set_river_inwater(
+    model::Model{N,L,V,R,W,T},
+    ssf_toriver,
+) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
+    @unpack lateral, vertical, network, config = model
     inds = network.index_river
-
-    @. lateral.river.inwater = (
-        ssf_toriver[inds] +
-        lateral.land.to_river[inds] +
-        # net_runoff_river
-        (
-            (
-                vertical.net_runoff_river[inds] *
-                network.land.xl[inds] *
-                network.land.yl[inds] *
-                0.001
-            ) / vertical.dt
+    do_water_demand = haskey(config.model, "water_demand")
+    if do_water_demand
+        @. lateral.river.inwater = (
+            ssf_toriver[inds] +
+            lateral.land.to_river[inds] +
+            # net_runoff_river
+            (vertical.net_runoff_river[inds] * network.land.area[inds] * 0.001) /
+            vertical.dt +
+            (lateral.river.allocation.nonirri_returnflow * 0.001 * network.river.area) /
+            vertical.dt
         )
-    )
+    else
+        @. lateral.river.inwater = (
+            ssf_toriver[inds] +
+            lateral.land.to_river[inds] +
+            # net_runoff_river
+            (vertical.net_runoff_river[inds] * network.land.area[inds] * 0.001) /
+            vertical.dt
+        )
+    end
 end
 
 """
@@ -1646,17 +1679,22 @@ Set `inwater` of the lateral land component for the `SbmGwfModel` type.
 """
 function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     @unpack lateral, vertical, network, config = model
-
+    do_water_demand = haskey(config.model, "water_demand")
     do_drains = get(config.model, "drains", false)::Bool
     drainflux = zeros(vertical.n)
     if do_drains
         drainflux[lateral.subsurface.drain.index] =
             -lateral.subsurface.drain.flux ./ tosecond(basetimestep)
     end
-
-    lateral.land.inwater .=
-        (vertical.net_runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.dt .+ drainflux
+    if do_water_demand
+        @. lateral.land.inwater =
+            (vertical.net_runoff + vertical.allocation.nonirri_returnflow) *
+            network.land.area *
+            0.001 / lateral.land.dt + drainflux
+    else
+        @. lateral.land.inwater =
+            (vertical.net_runoff * network.land.area * 0.001) / lateral.land.dt + drainflux
+    end
 end
 
 """
@@ -1665,10 +1703,17 @@ end
 Set `inwater` of the lateral land component for the `SbmModel` type.
 """
 function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
-    @unpack lateral, vertical, network = model
-    lateral.land.inwater .=
-        (vertical.net_runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.dt
+    @unpack lateral, vertical, network, config = model
+    do_water_demand = haskey(config.model, "water_demand")
+    if do_water_demand
+        @. lateral.land.inwater =
+            (vertical.net_runoff + vertical.allocation.nonirri_returnflow) *
+            network.land.area *
+            0.001 / lateral.land.dt
+    else
+        @. lateral.land.inwater =
+            (vertical.net_runoff * network.land.area * 0.001) / lateral.land.dt
+    end
 end
 
 """
@@ -1678,8 +1723,8 @@ Set `inwater` of the lateral land component, based on `runoff` of the `vertical`
 """
 function set_land_inwater(model)
     @unpack lateral, vertical, network = model
-    lateral.land.inwater .=
-        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.dt
+    @. lateral.land.inwater =
+        (vertical.runoff * network.land.area * 0.001) / lateral.land.dt
 end
 
 # Computation of inflow from the lateral components `land` and `subsurface` to water bodies
@@ -1804,13 +1849,10 @@ function surface_routing(
     @unpack lateral, vertical, network, clock = model
 
     @. lateral.land.runoff = (
-        (vertical.net_runoff / 1000.0) * (network.land.xl * network.land.yl) / vertical.dt +
+        (vertical.net_runoff / 1000.0) * network.land.area / vertical.dt +
         ssf_toriver +
         # net_runoff_river
-        (
-            (vertical.net_runoff_river * network.land.xl * network.land.yl * 0.001) /
-            vertical.dt
-        )
+        ((vertical.net_runoff_river * network.land.area * 0.001) / vertical.dt)
     )
     set_inflow_waterbody(model)
     update(lateral.land, lateral.river, network, julian_day(clock.time - clock.dt))
diff --git a/src/groundwater/aquifer.jl b/src/groundwater/aquifer.jl
index d650bf799..ea03659c7 100644
--- a/src/groundwater/aquifer.jl
+++ b/src/groundwater/aquifer.jl
@@ -95,6 +95,7 @@ transmissivity).
     specific_storage::Vector{T} | "m m-1 m-1" # [m m⁻¹ m⁻¹]
     storativity::Vector{T} | "m m-1" # [m m⁻¹]
     conductance::Vector{T} | "m2 d-1" # Confined aquifer conductance is constant
+    volume::Vector{T} | "m3" # total volume of water that can be released
 end
 
 
@@ -117,6 +118,7 @@ instead. Specific yield will vary roughly between 0.05 (clay) and 0.45 (peat)
     area::Vector{T} | "m2"
     specific_yield::Vector{T} | "m m-1" # [m m⁻¹]
     conductance::Vector{T} | "m2 d-1" #
+    volume::Vector{T} | "m3" # total volume of water that can be released
     f::Vector{T} | "-" # factor controlling the reduction of reference horizontal conductivity
     # Unconfined aquifer conductance is computed with degree of saturation (only when
     # conductivity_profile is set to "exponential")
@@ -158,6 +160,13 @@ function saturated_thickness(aquifer::ConfinedAquifer, index::Int)
     aquifer.top[index] - aquifer.bottom[index]
 end
 
+function saturated_thickness(aquifer::UnconfinedAquifer)
+    @. min(aquifer.top, aquifer.head) - aquifer.bottom
+end
+
+function saturated_thickness(aquifer::ConfinedAquifer)
+    @. aquifer.top - aquifer.bottom
+end
 
 """
     horizontal_conductance(i, j, nzi, aquifer, C)
@@ -355,6 +364,8 @@ function update(gwf, Q, dt, conductivity_profile)
     gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head
     # Make sure no heads ends up below an unconfined aquifer bottom
     gwf.aquifer.head .= minimum_head(gwf.aquifer)
+    gwf.aquifer.volume .=
+        saturated_thickness(gwf.aquifer) .* gwf.aquifer.area .* storativity(gwf.aquifer)
     return gwf
 end
 
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index aab5cbbd7..988dc4249 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -355,8 +355,7 @@ function initialize_hbv_model(config::Config)
         order = toposort,
         indices = inds,
         reverse_indices = rev_inds,
-        xl = xl,
-        yl = yl,
+        area = xl .* yl,
     )
     river = (
         graph = graph_riv,
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 0b6acf768..5060c842d 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -444,6 +444,19 @@ function initialize_storage(storfunc, area, waterlevel, sh)
     return storage
 end
 
+"Determine the water level depending on the storage function"
+function waterlevel(storfunc, area, storage, sh)
+    waterlevel = similar(area)
+    for i in eachindex(storage)
+        if storfunc[i] == 1
+            waterlevel[i] = storage[i] / area[i]
+        else
+            waterlevel[i] = interpolate_linear(storage[i], sh[i].S, sh[i].H)
+        end
+    end
+    return waterlevel
+end
+
 "Determine the maximum storage for lakes with a rating curve of type 1"
 function maximum_storage(storfunc, outflowfunc, area, sh, hq)
     maxstorage = Vector{Union{Float,Missing}}(missing, length(area))
diff --git a/src/sbm.jl b/src/sbm.jl
index 352ab5c08..60bfc6f24 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -35,6 +35,8 @@
     infiltcappath::Vector{T}
     # Soil infiltration capacity [mm Δt⁻¹]
     infiltcapsoil::Vector{T}
+    # Soil infiltration reduction factor (when soil is frozen) [-]
+    soilinfredu::Vector{T} | "-"
     # Maximum leakage [mm Δt⁻¹] from saturated zone
     maxleakage::Vector{T}
     # Fraction of open water (excluding rivers) [-]
@@ -43,6 +45,22 @@
     pathfrac::Vector{T} | "-"
     # Rooting depth [mm]
     rootingdepth::Vector{T} | "mm"
+    # Fraction of the root length density in each soil layer [-]
+    rootfraction::Vector{SVector{N,T}} | "-"
+    # Soil water pressure head h1 of the root water uptake reduction function (Feddes) [cm]
+    h1::Vector{T} | "cm"
+    # Soil water pressure head h2 of the root water uptake reduction function (Feddes) [cm]
+    h2::Vector{T} | "cm"
+    # Soil water pressure head h3_high of the root water uptake reduction function (Feddes) [cm]
+    h3_high::Vector{T} | "cm"
+    # Soil water pressure head h3_low of the root water uptake reduction function (Feddes) [cm]
+    h3_low::Vector{T} | "cm"
+    # Soil water pressure head h4 of the root water uptake reduction function (Feddes) [cm]
+    h4::Vector{T} | "cm"
+    # Calculated soil water pressure head h3 of the root water uptake reduction function (Feddes) [cm]
+    h3::Vector{T} | "cm"
+    # Root water uptake reduction at soil water pressure head h1 (0.0 or 1.0) [-]
+    alpha_h1::Vector{T} | "-"
     # Controls how roots are linked to water table [-]
     rootdistpar::Vector{T} | "-"
     # Parameter [mm] controlling capillary rise
@@ -206,6 +224,14 @@
     waterlevel_river::Vector{T} | "mm"
     # Total water storage (excluding floodplain volume, lakes and reservoirs) [mm]
     total_storage::Vector{T} | "mm"
+    # Water demand structs (of arrays)
+    paddy::Union{Paddy,Nothing} | "-" | 0
+    nonpaddy::Union{NonPaddy,Nothing} | "-" | 0
+    domestic::Union{NonIrrigationDemand,Nothing} | "-" | 0
+    livestock::Union{NonIrrigationDemand,Nothing} | "-" | 0
+    industry::Union{NonIrrigationDemand,Nothing} | "-" | 0
+    allocation::Union{AllocationLand,Nothing} | "-" | 0
+
 
     function SBM{T,N,M}(args...) where {T,N,M}
         equal_size_vectors(args)
@@ -348,33 +374,24 @@ function initialize_sbm(nc, config, riverfrac, inds)
         fill = 0.0,
     )
     # soil parameters
-    theta_s = ncread(
-        nc,
-        config,
-        "vertical.theta_s";
-        sel = inds,
-        defaults = 0.6,
-        type = Float,
-    )
-    theta_r = ncread(
-        nc,
-        config,
-        "vertical.theta_r";
-        sel = inds,
-        defaults = 0.01,
-        type = Float,
-    )
+    theta_s =
+        ncread(nc, config, "vertical.theta_s"; sel = inds, defaults = 0.6, type = Float)
+    theta_r =
+        ncread(nc, config, "vertical.theta_r"; sel = inds, defaults = 0.01, type = Float)
     kv_0 =
-        ncread(
-            nc,
-            config,
-            "vertical.kv_0";
-            sel = inds,
-            defaults = 3000.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(nc, config, "vertical.kv_0"; sel = inds, defaults = 3000.0, type = Float) .*
+        (dt / basetimestep)
     f = ncread(nc, config, "vertical.f"; sel = inds, defaults = 0.001, type = Float)
-    hb = ncread(nc, config, "vertical.hb"; sel = inds, defaults = 10.0, type = Float)
+    hb = ncread(nc, config, "vertical.hb"; sel = inds, defaults = -10.0, type = Float)
+    h1 = ncread(nc, config, "vertical.h1"; sel = inds, defaults = 0.0, type = Float)
+    h2 = ncread(nc, config, "vertical.h2"; sel = inds, defaults = -100.0, type = Float)
+    h3_high =
+        ncread(nc, config, "vertical.h3_high"; sel = inds, defaults = -400.0, type = Float)
+    h3_low =
+        ncread(nc, config, "vertical.h3_low"; sel = inds, defaults = -1000.0, type = Float)
+    h4 = ncread(nc, config, "vertical.h4"; sel = inds, defaults = -15849.0, type = Float)
+    alpha_h1 =
+        ncread(nc, config, "vertical.alpha_h1"; sel = inds, defaults = 1.0, type = Float)
     soilthickness = ncread(
         nc,
         config,
@@ -455,6 +472,8 @@ function initialize_sbm(nc, config, riverfrac, inds)
         defaults = 750.0,
         type = Float,
     )
+    # correct rooting depth for soilthickness
+    rootingdepth = @. min(0.99 * soilthickness, rootingdepth)
     rootdistpar = ncread(
         nc,
         config,
@@ -512,6 +531,42 @@ function initialize_sbm(nc, config, riverfrac, inds)
             s_layers[:, i] = pushfirst(cumsum(SVector(soilthickness[i])), 0.0)
         end
     end
+
+    if length(config_thicknesslayers) > 0
+        # root fraction read from nc file, in case of multiple soil layers and TOML file
+        # includes "vertical.rootfraction"
+        if haskey(config.input.vertical, "rootfraction")
+            rootfraction = ncread(
+                nc,
+                config,
+                "vertical.rootfraction";
+                sel = inds,
+                optional = false,
+                type = Float,
+                dimname = :layer,
+            )
+        else
+            # default root fraction in case of multiple soil layers
+            rootfraction = zeros(Float, maxlayers, n)
+            for i = 1:n
+                if rootingdepth[i] > 0.0
+                    for k = 1:maxlayers
+                        if (rootingdepth[i] - s_layers[k, i]) >= act_thickl[k, i]
+                            rootfraction[k, i] = act_thickl[k, i] / rootingdepth[i]
+                        else
+                            rootfraction[k, i] =
+                                max(rootingdepth[i] - s_layers[k, i], 0.0) / rootingdepth[i]
+                        end
+                    end
+                end
+            end
+        end
+    else
+        # for the case of 1 soil layer
+        rootfraction = ones(Float, maxlayers, n)
+    end
+
+    # needed for derived parameters below
     act_thickl = svectorscopy(act_thickl, Val{maxlayers}())
 
     # copied to array of sarray below
@@ -575,7 +630,15 @@ function initialize_sbm(nc, config, riverfrac, inds)
               """)
     end
 
-    sbm = SBM{Float,maxlayers,maxlayers + 1}(
+    # water demand and irrigation options
+    do_water_demand = haskey(config.model, "water_demand")
+    domestic = do_water_demand ? get(config.model.water_demand, "domestic", false) : false
+    industry = do_water_demand ? get(config.model.water_demand, "industry", false) : false
+    livestock = do_water_demand ? get(config.model.water_demand, "livestock", false) : false
+    paddy = do_water_demand ? get(config.model.water_demand, "paddy", false) : false
+    nonpaddy = do_water_demand ? get(config.model.water_demand, "nonpaddy", false) : false
+
+    sbm = SBM(
         dt = tosecond(dt),
         maxlayers = maxlayers,
         n = n,
@@ -589,15 +652,24 @@ function initialize_sbm(nc, config, riverfrac, inds)
         kv = svectorscopy(kv, Val{maxlayers}()),
         kvfrac = svectorscopy(kvfrac, Val{maxlayers}()),
         hb = hb,
+        h1 = h1,
+        h2 = h2,
+        h3_high = h3_high,
+        h3_low = h3_low,
+        h4 = h4,
+        h3 = fill(mv, n),
+        alpha_h1 = alpha_h1,
         soilthickness = soilthickness,
         act_thickl = act_thickl,
         sumlayers = sumlayers,
         infiltcappath = infiltcappath,
         infiltcapsoil = infiltcapsoil,
+        soilinfredu = fill(Float(1), n),
         maxleakage = maxleakage,
         waterfrac = max.(waterfrac .- riverfrac, Float(0.0)),
         pathfrac = pathfrac,
         rootingdepth = rootingdepth,
+        rootfraction = svectorscopy(rootfraction, Val{maxlayers}()),
         rootdistpar = rootdistpar,
         cap_hmax = cap_hmax,
         cap_n = cap_n,
@@ -678,9 +750,18 @@ function initialize_sbm(nc, config, riverfrac, inds)
         swood = swood,
         kext = kext,
         leaf_area_index = fill(mv, n),
+        # water level land and river domain
         waterlevel_land = fill(mv, n),
         waterlevel_river = zeros(Float, n), #set to zero to account for cells outside river domain
         total_storage = zeros(Float, n), # Set the total water storage from initialized values
+        # water demand
+        paddy = paddy ? initialize_paddy(nc, config, inds, dt) : nothing,
+        nonpaddy = nonpaddy ? initialize_nonpaddy(nc, config, inds, dt) : nothing,
+        domestic = domestic ? initialize_domestic_demand(nc, config, inds, dt) : nothing,
+        industry = industry ? initialize_industry_demand(nc, config, inds, dt) : nothing,
+        livestock = livestock ? initialize_livestock_demand(nc, config, inds, dt) : nothing,
+        allocation = do_water_demand ? initialize_allocation_land(nc, config, inds) :
+                     nothing,
     )
 
     return sbm
@@ -747,6 +828,9 @@ function update_until_snow(sbm::SBM, config)
                 sbm.whc[i],
             )
         end
+
+        h3 = feddes_h3(sbm.h3_high[i], sbm.h3_low[i], pottrans, Second(sbm.dt))
+
         # update the outputs and states
         sbm.e_r[i] = e_r
         sbm.cmax[i] = cmax
@@ -756,6 +840,7 @@ function update_until_snow(sbm::SBM, config)
         sbm.stemflow[i] = stemflow
         sbm.throughfall[i] = throughfall
         sbm.pottrans[i] = pottrans
+        sbm.h3[i] = h3
         if modelsnow
             sbm.snow[i] = snow
             sbm.snowwater[i] = snowwater
@@ -807,10 +892,11 @@ function update_until_recharge(sbm::SBM, config)
 
         runoff_river = min(1.0, sbm.riverfrac[i]) * avail_forinfilt
         runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt
+        if !isnothing(sbm.paddy) || !isnothing(sbm.nonpaddy)
+            avail_forinfilt = avail_forinfilt + sbm.allocation.irri_alloc[i]
+        end
         avail_forinfilt = max(avail_forinfilt - runoff_river - runoff_land, 0.0)
 
-        rootingdepth = min(sbm.soilthickness[i] * 0.99, sbm.rootingdepth[i])
-
         ae_openw_r = min(
             sbm.waterlevel_river[i] * sbm.riverfrac[i],
             sbm.riverfrac[i] * sbm.potential_evaporation[i],
@@ -828,22 +914,36 @@ function update_until_recharge(sbm::SBM, config)
         )
         potsoilevap = soilevap_fraction * sbm.potential_evaporation[i]
 
+        if !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i]
+            evap_paddy_water = min(sbm.paddy.h[i], potsoilevap)
+            sbm.paddy.h[i] -= evap_paddy_water
+            potsoilevap -= evap_paddy_water
+            avail_forinfilt += sbm.paddy.h[i] # allow infiltration of paddy water
+        else
+            evap_paddy_water = 0.0
+        end
+
         # Calculate the initial capacity of the unsaturated store
         ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
 
         # Calculate the infiltration flux into the soil column
-        infiltsoilpath, infiltsoil, infiltpath, soilinf, pathinf, infiltexcess =
-            infiltration(
-                avail_forinfilt,
-                sbm.pathfrac[i],
-                sbm.cf_soil[i],
-                sbm.tsoil[i],
-                sbm.infiltcapsoil[i],
-                sbm.infiltcappath[i],
-                ustorecapacity,
-                modelsnow,
-                soilinfreduction,
-            )
+        infiltsoilpath,
+        infiltsoil,
+        infiltpath,
+        soilinf,
+        pathinf,
+        infiltexcess,
+        soilinfredu = infiltration(
+            avail_forinfilt,
+            sbm.pathfrac[i],
+            sbm.cf_soil[i],
+            sbm.tsoil[i],
+            sbm.infiltcapsoil[i],
+            sbm.infiltcappath[i],
+            ustorecapacity,
+            modelsnow,
+            soilinfreduction,
+        )
 
 
         usl, n_usl = set_layerthickness(sbm.zi[i], sbm.sumlayers[i], sbm.act_thickl[i])
@@ -933,31 +1033,73 @@ function update_until_recharge(sbm::SBM, config)
         soilevap = soilevapunsat + soilevapsat
         satwaterdepth = sbm.satwaterdepth[i] - soilevapsat
 
-        # transpiration from saturated store
-        wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
-        actevapsat = min(sbm.pottrans[i] * wetroots, satwaterdepth)
-        satwaterdepth = satwaterdepth - actevapsat
-        restpottrans = sbm.pottrans[i] - actevapsat
-
-        # actual transpiration from ustore
+        # actual transpiration, first from ustore, potential transpiration is partitioned over
+        # depth based on the rootfraction
         actevapustore = 0.0
+        rootfraction_unsat = 0.0
         for k = 1:n_usl
-            ustorelayerdepth, actevapustore, restpottrans = acttransp_unsat_sbm(
-                rootingdepth,
-                usld[k],
-                sbm.sumlayers[i][k],
-                restpottrans,
-                actevapustore,
-                sbm.c[i][k],
-                usl[k],
+            vwc = max(usld[k] / usl[k], Float(0.0000001))
+            head = head_brooks_corey(
+                vwc,
                 sbm.theta_s[i],
                 sbm.theta_r[i],
+                sbm.c[i][k],
                 sbm.hb[i],
-                ust,
             )
+            alpha = rwu_reduction_feddes(
+                head,
+                sbm.h1[i],
+                sbm.h2[i],
+                sbm.h3[i],
+                sbm.h4[i],
+                sbm.alpha_h1[i],
+            )
+            # availcap is fraction of soil layer containing roots
+            # if `ust` is `true`, the whole unsaturated store is available for transpiration
+            if ust
+                availcap = usld[k] * 0.99
+            else
+                availcap =
+                    min(1.0, max(0.0, (sbm.rootingdepth[i] - sbm.sumlayers[i][k]) / usl[k]))
+            end
+            maxextr = usld[k] * availcap
+            # the rootfraction is valid for the root length in a soil layer, if zi decreases the root length
+            # the rootfraction needs to be adapted
+            if k == n_usl && sbm.zi[i] < sbm.rootingdepth[i]
+                rootlength =
+                    min(sbm.act_thickl[i][k], sbm.rootingdepth[i] - sbm.sumlayers[i][k])
+                rootfraction_act = sbm.rootfraction[i][k] * (usl[k] / rootlength)
+            else
+                rootfraction_act = sbm.rootfraction[i][k]
+            end
+            actevapustore_layer = min(alpha * rootfraction_act * sbm.pottrans[i], maxextr)
+            rootfraction_unsat = rootfraction_unsat + rootfraction_act
+            ustorelayerdepth = usld[k] - actevapustore_layer
+            actevapustore = actevapustore + actevapustore_layer
             usld = setindex(usld, ustorelayerdepth, k)
         end
 
+        # transpiration from saturated store
+        wetroots = scurve(sbm.zi[i], sbm.rootingdepth[i], Float(1.0), sbm.rootdistpar[i])
+        alpha = rwu_reduction_feddes(
+            Float(0.0),
+            sbm.h1[i],
+            sbm.h2[i],
+            sbm.h3[i],
+            sbm.h4[i],
+            sbm.alpha_h1[i],
+        )
+        # include remaining root fraction if rooting depth is below water table zi
+        if sbm.zi[i] >= sbm.rootingdepth[i]
+            f_roots = wetroots
+            restevap = sbm.pottrans[i] - actevapustore
+        else
+            f_roots = wetroots * (1.0 - rootfraction_unsat)
+            restevap = sbm.pottrans[i]
+        end
+        actevapsat = min(restevap * f_roots * alpha, satwaterdepth)
+        satwaterdepth = satwaterdepth - actevapsat
+
         # check soil moisture balance per layer
         du = 0.0
         for k = n_usl:-1:1
@@ -990,7 +1132,7 @@ function update_until_recharge(sbm::SBM, config)
                 sbm.soilwatercapacity[i] - satwaterdepth - sum(@view usld[1:sbm.nlayers[i]])
             maxcapflux = max(0.0, min(ksat, actevapustore, ustorecapacity, satwaterdepth))
 
-            if sbm.zi[i] > rootingdepth
+            if sbm.zi[i] > sbm.rootingdepth[i]
                 capflux =
                     maxcapflux * pow(
                         1.0 - min(sbm.zi[i], sbm.cap_hmax[i]) / (sbm.cap_hmax[i]),
@@ -1002,8 +1144,10 @@ function update_until_recharge(sbm::SBM, config)
 
             netcapflux = capflux
             for k = n_usl:-1:1
-                toadd =
-                    min(netcapflux, max(usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k], 0.0))
+                toadd = min(
+                    netcapflux,
+                    max(usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k], 0.0),
+                )
                 usld = setindex(usld, usld[k] + toadd, k)
                 netcapflux = netcapflux - toadd
                 actcapflux = actcapflux + toadd
@@ -1022,7 +1166,13 @@ function update_until_recharge(sbm::SBM, config)
         # recharge (mm) for saturated zone
         recharge = (transfer - actcapflux - actleakage - actevapsat - soilevapsat)
         transpiration = actevapsat + actevapustore
-        actevap = soilevap + transpiration + ae_openw_r + ae_openw_l + sbm.interception[i]
+        actevap =
+            soilevap +
+            transpiration +
+            ae_openw_r +
+            ae_openw_l +
+            sbm.interception[i] +
+            evap_paddy_water
 
         # update the outputs and states
         sbm.n_unsatlayers[i] = n_usl
@@ -1053,6 +1203,7 @@ function update_until_recharge(sbm::SBM, config)
         sbm.rainfallplusmelt[i] = rainfallplusmelt
         sbm.infiltsoilpath[i] = infiltsoilpath
         sbm.satwaterdepth[i] = satwaterdepth
+        sbm.soilinfredu[i] = soilinfredu
         if modelsnow
             if modelglacier
                 sbm.snow[i] = snow
@@ -1083,12 +1234,23 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
 
         ustoredepth = sum(@view usld[1:n_usl])
 
-        runoff =
-            exfiltustore +
-            exfiltsatwater[i] +
-            sbm.excesswater[i] +
-            sbm.runoff_land[i] +
-            sbm.infiltexcess[i]
+        if !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i]
+            paddy_h_add =
+                exfiltustore +
+                exfiltsatwater[i] +
+                sbm.excesswater[i] +
+                sbm.runoff_land[i] +
+                sbm.infiltexcess[i]
+            runoff = max(paddy_h_add - sbm.paddy.h_max[i], 0.0)
+            sbm.paddy.h[i] = paddy_h_add - runoff
+        else
+            runoff =
+                exfiltustore +
+                exfiltsatwater[i] +
+                sbm.excesswater[i] +
+                sbm.runoff_land[i] +
+                sbm.infiltexcess[i]
+        end
 
         # volumetric water content per soil layer and root zone
         vwc = sbm.vwc[i]
@@ -1097,7 +1259,10 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
             if k <= n_usl
                 vwc = setindex(
                     vwc,
-                    (usld[k] + (sbm.act_thickl[i][k] - usl[k]) * (sbm.theta_s[i] - sbm.theta_r[i])) / sbm.act_thickl[i][k] + sbm.theta_r[i],
+                    (
+                        usld[k] +
+                        (sbm.act_thickl[i][k] - usl[k]) * (sbm.theta_s[i] - sbm.theta_r[i])
+                    ) / sbm.act_thickl[i][k] + sbm.theta_r[i],
                     k,
                 )
             else
@@ -1114,7 +1279,8 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
                 usld[k]
         end
 
-        rootstore_sat = max(0.0, sbm.rootingdepth[i] - zi[i]) * (sbm.theta_s[i] - sbm.theta_r[i])
+        rootstore_sat =
+            max(0.0, sbm.rootingdepth[i] - zi[i]) * (sbm.theta_s[i] - sbm.theta_r[i])
         rootstore = rootstore_sat + rootstore_unsat
         vwc_root = rootstore / sbm.rootingdepth[i] + sbm.theta_r[i]
         vwc_percroot = (vwc_root / sbm.theta_s[i]) * 100.0
@@ -1147,10 +1313,8 @@ Takes the following parameters:
     The vertical concept (SBM struct)
 - river_network:
     The indices of the river cells in relation to the active cells, i.e. model.network.index_river
-- cell_xsize:
-    Size in X direction of the cells acquired from model.network.land.xl
-- cell_ysize:
-    Size in Y direction of the cells acquired from model.network.land.yl
+- area:
+    Area of the cells acquired from model.network.land.area
 - river_routing:
     The river routing struct, i.e. model.lateral.river
 - land_routing:
@@ -1159,8 +1323,7 @@ Takes the following parameters:
 function update_total_water_storage(
     sbm::SBM,
     river_network,
-    cell_xsize,
-    cell_ysize,
+    area,
     river_routing,
     land_routing,
 )
@@ -1175,19 +1338,23 @@ function update_total_water_storage(
         (
             river_routing.h_av[1:nriv] .* river_routing.width[1:nriv] .*
             river_routing.dl[1:nriv]
-        ) ./ (cell_xsize[river_network] .* cell_ysize[river_network]) * 1000 # Convert to mm
+        ) ./ (area[river_network]) * 1000 # Convert to mm
     )
 
     # Chunk the data for parallel computing
     threaded_foreach(1:sbm.n, basesize = 1000) do i
 
+        # Paddy water depth
+        paddy_h = isnothing(sbm.paddy) ? 0.0 : sbm.paddy.h[i]
+
         # Cumulate per vertical type
         # Maybe re-categorize in the future
         surface = (
             sbm.glacierstore[i] * sbm.glacierfrac[i] +
             sbm.snow[i] +
             sbm.snowwater[i] +
-            sbm.canopystorage[i]
+            sbm.canopystorage[i] +
+            paddy_h
         )
         sub_surface = sbm.ustoredepth[i] + sbm.satwaterdepth[i]
         lateral = (
@@ -1198,3 +1365,104 @@ function update_total_water_storage(
         sbm.total_storage[i] += (surface + sub_surface + lateral)
     end
 end
+
+"""
+    update_water_demand(sbm::SBM)
+
+Update water demand for vertical `SBM` concept for a single timestep. Water demand is
+computed for sectors `industry`, `domestic` and `livestock`, and `paddy` rice fields and
+`nonpaddy` (other crop) fields.
+
+Gross water demand for irrigation `irri_demand_gross` and non-irrigation
+`nonirri_demand_gross`, and total gross water demand `total_gross_demand` are updated as
+part of `SBM` water allocation (`allocation`)
+"""
+function update_water_demand(sbm::SBM)
+    for i = 1:sbm.n
+
+        industry_dem = update_non_irrigation_demand(sbm.industry, i)
+        domestic_dem = update_non_irrigation_demand(sbm.domestic, i)
+        livestock_dem = update_non_irrigation_demand(sbm.livestock, i)
+
+        irri_dem_gross = 0.0
+        if !isnothing(sbm.nonpaddy) && sbm.nonpaddy.irrigation_areas[i]
+            if sbm.nonpaddy.irrigation_trigger[i]
+                usl, _ = set_layerthickness(sbm.zi[i], sbm.sumlayers[i], sbm.act_thickl[i])
+                for k = 1:sbm.n_unsatlayers[i]
+                    # compute water demand only for root zone through root fraction per layer
+                    rootfrac = min(
+                        1.0,
+                        (max(0.0, sbm.rootingdepth[i] - sbm.sumlayers[i][k]) / usl[k]),
+                    )
+                    # vwc_f and vwc_h3 can be precalculated.
+                    vwc_fc = vwc_brooks_corey(
+                        -100.0,
+                        sbm.hb[i],
+                        sbm.theta_s[i],
+                        sbm.theta_r[i],
+                        sbm.c[i][k],
+                    )
+                    vwc_h3 = vwc_brooks_corey(
+                        sbm.h3[i],
+                        sbm.hb[i],
+                        sbm.theta_s[i],
+                        sbm.theta_r[i],
+                        sbm.c[i][k],
+                    )
+                    depletion =
+                        (vwc_fc * usl[k]) -
+                        (sbm.ustorelayerdepth[i][k] + sbm.theta_r[i] * usl[k])
+                    depletion *= rootfrac
+                    raw = (vwc_fc - vwc_h3) * usl[k] # readily available water
+                    raw *= rootfrac
+
+                    # check if maximum irrigation rate has been applied at the previous time step.
+                    max_irri_rate_applied =
+                        sbm.nonpaddy.demand_gross[i] ==
+                        sbm.nonpaddy.maximum_irrigation_rate[i]
+                    if depletion >= raw # start irrigation
+                        irri_dem_gross += depletion
+                        # add depletion to irrigation gross demand when the maximum irrigation rate has been 
+                        # applied at the previous time step (to get volumetric water content at field capacity)
+                    elseif depletion > 0.0 && max_irri_rate_applied # continue irrigation
+                        irri_dem_gross += depletion
+                    end
+                end
+                # limit irrigation demand to infiltration capacity    
+                infiltration_capacity =
+                    sbm.soilinfredu[i] * (1.0 - sbm.pathfrac[i]) * sbm.infiltcapsoil[i]
+                irri_dem_gross = min(irri_dem_gross, infiltration_capacity)
+                irri_dem_gross /= sbm.nonpaddy.irrigation_efficiency[i]
+                # limit irrigation demand to the maximum irrigation rate
+                irri_dem_gross =
+                    min(irri_dem_gross, sbm.nonpaddy.maximum_irrigation_rate[i])
+            else
+                irri_dem_gross = 0.0
+            end
+            sbm.nonpaddy.demand_gross[i] = irri_dem_gross
+        elseif !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i]
+            if sbm.paddy.irrigation_trigger[i]
+                # check if maximum irrigation rate has been applied at the previous time step.
+                max_irri_rate_applied =
+                    sbm.paddy.demand_gross[i] == sbm.paddy.maximum_irrigation_rate[i]
+                # start irrigation
+                if sbm.paddy.h[i] < sbm.paddy.h_min[i]
+                    irr_depth_paddy = sbm.paddy.h_opt[i] - sbm.paddy.h[i]
+                elseif sbm.paddy.h[i] < sbm.paddy.h_opt[i] && max_irri_rate_applied # continue irrigation
+                    irr_depth_paddy = sbm.paddy.h_opt[i] - sbm.paddy.h[i]
+                else
+                    irr_depth_paddy = 0.0
+                end
+                irri_dem_gross += irr_depth_paddy / sbm.paddy.irrigation_efficiency[i]
+                # limit irrigation demand to the maximum irrigation rate
+                irri_dem_gross = min(irri_dem_gross, sbm.paddy.maximum_irrigation_rate[i])
+            end
+            sbm.paddy.demand_gross[i] = irri_dem_gross
+        end
+        # update gross water demands 
+        sbm.allocation.irri_demand_gross[i] = irri_dem_gross
+        sbm.allocation.nonirri_demand_gross[i] = industry_dem + domestic_dem + livestock_dem
+        sbm.allocation.total_gross_demand[i] =
+            irri_dem_gross + industry_dem + domestic_dem + livestock_dem
+    end
+end
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 6a76b5bdf..faf861171 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -41,6 +41,7 @@ function initialize_sbm_gwf_model(config::Config)
     )::String
     land_routing =
         get_options(config.model, "land_routing", routing_options, "kinematic-wave")::String
+    do_water_demand = haskey(config.model, "water_demand")
 
     nc = NCDataset(static_path)
 
@@ -71,12 +72,12 @@ function initialize_sbm_gwf_model(config::Config)
     xl, yl = cell_lengths(y, cellength, sizeinmetres)
     riverfrac = river_fraction(river, riverlength, riverwidth, xl, yl)
 
-    # initialize vertical SBM concept
-    sbm = initialize_sbm(nc, config, riverfrac, inds)
-
     inds_riv, rev_inds_riv = active_indices(river_2d, 0)
     nriv = length(inds_riv)
 
+    # initialize vertical SBM concept
+    sbm = initialize_sbm(nc, config, riverfrac, inds)
+
     # reservoirs
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
@@ -118,6 +119,18 @@ function initialize_sbm_gwf_model(config::Config)
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
     frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
+    inds_allocation_areas = Vector{Int}[]
+    inds_riv_allocation_areas = Vector{Int}[]
+    if do_water_demand
+        areas = unique(sbm.allocation.areas)
+        for a in areas
+            area_index = findall(x -> x == a, sbm.allocation.areas)
+            push!(inds_allocation_areas, area_index)
+            area_riv_index = findall(x -> x == a, sbm.allocation.areas[index_river])
+            push!(inds_riv_allocation_areas, area_riv_index)
+        end
+    end
+
     if land_routing == "kinematic-wave"
         olf = initialize_surfaceflow_land(
             nc,
@@ -227,28 +240,29 @@ function initialize_sbm_gwf_model(config::Config)
         get(config.input.lateral.subsurface, "conductivity_profile", "uniform")::String
 
     connectivity = Connectivity(inds, rev_inds, xl, yl)
-    initial_head = altitude .- Float(0.10) # cold state for groundwater head
+    initial_head = altitude .- sbm.zi / 1000.0 # cold state for groundwater head based on SBM zi
     initial_head[index_river] = altitude[index_river]
 
     if do_constanthead
         initial_head[constant_head.index] = constant_head.head
     end
 
+    bottom = altitude .- sbm.soilthickness ./ Float(1000.0)
+    area = xl .* yl
+    volume = @. (min(altitude, initial_head) - bottom) * area * specific_yield # total volume than can be released
+
     aquifer = UnconfinedAquifer(
         initial_head,
         conductivity,
         altitude,
-        altitude .- sbm.soilthickness ./ Float(1000.0),
-        xl .* yl,
+        bottom,
+        area,
         specific_yield,
         zeros(Float, connectivity.nconnection),  # conductance
+        volume,
         gwf_f,
     )
 
-    # reset zi and satwaterdepth with groundwater head from unconfined aquifer
-    sbm.zi .= (altitude .- initial_head) .* 1000.0
-    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* (sbm.theta_s .- sbm.theta_r)
-
     # river boundary of unconfined aquifer
     infiltration_conductance = ncread(
         nc,
@@ -419,10 +433,10 @@ function initialize_sbm_gwf_model(config::Config)
             order = toposort,
             indices = inds,
             reverse_indices = rev_inds,
-            xl = xl,
-            yl = yl,
+            area = xl .* yl,
             slope = landslope,
             altitude = altitude,
+            indices_allocation_areas = inds_allocation_areas,
         )
     elseif land_routing == "local-inertial"
         land = (
@@ -430,34 +444,66 @@ function initialize_sbm_gwf_model(config::Config)
             order = toposort,
             indices = inds,
             reverse_indices = rev_inds,
-            xl = xl,
-            yl = yl,
+            area = xl .* yl,
             slope = landslope,
             altitude = altitude,
             index_river = index_river_nf,
             staggered_indices = indices,
+            indices_allocation_areas = inds_allocation_areas,
         )
     end
+    if do_water_demand
+        # exclude waterbodies for local surface and ground water abstraction
+        inds_riv_2d = copy(rev_inds_riv)
+        inds_2d = ones(Bool, modelsize_2d)
+        if !isempty(reservoir)
+            inds_cov = collect(Iterators.flatten(reservoir.indices_coverage))
+            inds_riv_2d[inds_cov] .= 0
+            inds_2d[inds_cov] .= 0
+        end
+        if !isempty(lake)
+            inds_cov = collect(Iterators.flatten(lake.indices_coverage))
+            inds_riv_2d[inds_cov] .= 0
+            inds_2d[inds_cov] .= 0
+        end
+        land = merge(land, (index_river_wb = inds_riv_2d[inds], index_wb = inds_2d[inds]))
+    end
     if river_routing == "kinematic-wave"
         river = (
             graph = graph_riv,
             indices = inds_riv,
             reverse_indices = rev_inds_riv,
+            # reservoir and lake index
+            reservoir_index = resindex,
+            lake_index = lakeindex,
+            reservoir_index_f = filter(x -> x ≠ 0, resindex),
+            lake_index_f = filter(x -> x ≠ 0, lakeindex),
             # specific for kinematic_wave
             upstream_nodes = filter_upsteam_nodes(graph_riv, pits[inds_riv]),
             subdomain_order = subriv_order,
             topo_subdomain = topo_subriv,
             indices_subdomain = indices_subriv,
             order = toposort_riv,
+            # water allocation areas
+            indices_allocation_areas = inds_riv_allocation_areas,
+            area = xl[index_river] .* yl[index_river],
         )
     elseif river_routing == "local-inertial"
         river = (
             graph = graph_riv,
             indices = inds_riv,
             reverse_indices = rev_inds_riv,
+            # reservoir and lake index
+            reservoir_index = resindex,
+            lake_index = lakeindex,
+            reservoir_index_f = filter(x -> x ≠ 0, resindex),
+            lake_index_f = filter(x -> x ≠ 0, lakeindex),
             # specific for local-inertial
             nodes_at_link = nodes_at_link,
             links_at_node = adjacent_links_at_node(graph_riv, nodes_at_link),
+            # water allocation areas
+            indices_allocation_areas = inds_riv_allocation_areas,
+            area = xl[index_river] .* yl[index_river],
         )
     end
 
@@ -482,6 +528,7 @@ end
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     @unpack lateral, vertical, network, clock, config = model
 
+    do_water_demand = haskey(config.model, "water_demand")
     inds_riv = network.index_river
 
     # extract water levels h_av [m] from the land and river domains
@@ -503,6 +550,12 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
         )
     end
 
+    # optional water demand and allocation
+    if do_water_demand
+        update_water_demand(vertical)
+        update_water_allocation(model)
+    end
+
     # update vertical sbm concept until recharge [mm]
     update_until_recharge(vertical, config)
 
@@ -523,7 +576,11 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     Q = zeros(vertical.n)
     # exchange of recharge between vertical sbm concept and groundwater flow domain
     # recharge rate groundwater is required in units [m d⁻¹]
-    lateral.subsurface.recharge.rate .= vertical.recharge ./ 1000.0 .* (1.0 / dt_sbm)
+    @. lateral.subsurface.recharge.rate = vertical.recharge / 1000.0 * (1.0 / dt_sbm)
+    if do_water_demand
+        @. lateral.subsurface.recharge.rate -=
+            vertical.allocation.act_groundwater_abst / 1000.0 * (1.0 / dt_sbm)
+    end
     # update groundwater domain
     update(lateral.subsurface.flow, Q, dt_sbm, conductivity_profile)
 
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index bc238bab9..ec1032894 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -33,6 +33,7 @@ function initialize_sbm_model(config::Config)
     )::String
     land_routing =
         get_options(config.model, "land_routing", routing_options, "kinematic-wave")::String
+    do_water_demand = haskey(config.model, "water_demand")
 
     snow = get(config.model, "snow", false)::Bool
     reservoirs = do_reservoirs
@@ -69,11 +70,11 @@ function initialize_sbm_model(config::Config)
     xl, yl = cell_lengths(y, cellength, sizeinmetres)
     riverfrac = river_fraction(river, riverlength, riverwidth, xl, yl)
 
-    sbm = initialize_sbm(nc, config, riverfrac, inds)
-
     inds_riv, rev_inds_riv = active_indices(river_2d, 0)
     nriv = length(inds_riv)
 
+    sbm = initialize_sbm(nc, config, riverfrac, inds)
+
     # reservoirs
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
@@ -149,6 +150,7 @@ function initialize_sbm_model(config::Config)
             ssfin = fill(mv, n),
             ssfmax = fill(mv, n),
             to_river = zeros(n),
+            volume = (sbm.theta_s .- sbm.theta_r) .* (soilthickness .- zi) .* (xl .* yl),
         )
         # update variables `ssf`, `ssfmax` and `kh` (layered profile) based on ksat_profile
         ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
@@ -182,6 +184,18 @@ function initialize_sbm_model(config::Config)
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
     frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
+    inds_allocation_areas = Vector{Int}[]
+    inds_riv_allocation_areas = Vector{Int}[]
+    if do_water_demand
+        areas = unique(sbm.allocation.areas)
+        for a in areas
+            area_index = findall(x -> x == a, sbm.allocation.areas)
+            push!(inds_allocation_areas, area_index)
+            area_riv_index = findall(x -> x == a, sbm.allocation.areas[index_river])
+            push!(inds_riv_allocation_areas, area_riv_index)
+        end
+    end
+
     if land_routing == "kinematic-wave"
         olf = initialize_surfaceflow_land(
             nc,
@@ -338,33 +352,65 @@ function initialize_sbm_model(config::Config)
         order = toposort,
         indices = inds,
         reverse_indices = rev_inds,
-        xl,
-        yl,
+        area = xl .* yl,
         slope = landslope,
+        indices_allocation_areas = inds_allocation_areas,
     )
     if land_routing == "local-inertial"
         land = merge(land, (index_river = index_river_nf, staggered_indices = indices))
     end
+    if do_water_demand
+        # exclude waterbodies for local surface and ground water abstraction
+        inds_riv_2d = copy(rev_inds_riv)
+        inds_2d = ones(Bool, modelsize_2d)
+        if !isempty(reservoir)
+            inds_cov = collect(Iterators.flatten(reservoir.indices_coverage))
+            inds_riv_2d[inds_cov] .= 0
+            inds_2d[inds_cov] .= 0
+        end
+        if !isempty(lake)
+            inds_cov = collect(Iterators.flatten(lake.indices_coverage))
+            inds_riv_2d[inds_cov] .= 0
+            inds_2d[inds_cov] .= 0
+        end
+        land = merge(land, (index_river_wb = inds_riv_2d[inds], index_wb = inds_2d[inds]))
+    end
     if river_routing == "kinematic-wave"
         river = (
             graph = graph_riv,
             indices = inds_riv,
             reverse_indices = rev_inds_riv,
+            # reservoir and lake index
+            reservoir_index = resindex,
+            lake_index = lakeindex,
+            reservoir_index_f = filter(x -> x ≠ 0, resindex),
+            lake_index_f = filter(x -> x ≠ 0, lakeindex),
             # specific for kinematic_wave
             upstream_nodes = filter_upsteam_nodes(graph_riv, pits[inds_riv]),
             subdomain_order = subriv_order,
             topo_subdomain = topo_subriv,
             indices_subdomain = indices_subriv,
             order = toposort_riv,
+            # water allocation areas
+            indices_allocation_areas = inds_riv_allocation_areas,
+            area = xl[index_river] .* yl[index_river],
         )
     elseif river_routing == "local-inertial"
         river = (
             graph = graph_riv,
             indices = inds_riv,
             reverse_indices = rev_inds_riv,
+            # reservoir and lake index
+            reservoir_index = resindex,
+            lake_index = lakeindex,
+            reservoir_index_f = filter(x -> x ≠ 0, resindex),
+            lake_index_f = filter(x -> x ≠ 0, lakeindex),
             # specific for local-inertial
             nodes_at_link = nodes_at_link,
             links_at_node = adjacent_links_at_node(graph_riv, nodes_at_link),
+            # water allocation areas
+            indices_allocation_areas = inds_riv_allocation_areas,
+            area = xl[index_river] .* yl[index_river],
         )
     end
 
@@ -389,11 +435,15 @@ end
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
 
     @unpack lateral, vertical, network, clock, config = model
+    do_water_demand = haskey(config.model, "water_demand")
     ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
 
     model = update_until_recharge(model)
     # exchange of recharge between vertical sbm concept and subsurface flow domain
     lateral.subsurface.recharge .= vertical.recharge ./ 1000.0
+    if do_water_demand
+        @. lateral.subsurface.recharge -= vertical.allocation.act_groundwater_abst / 1000.0
+    end
     lateral.subsurface.recharge .*= lateral.subsurface.dw
     lateral.subsurface.zi .= vertical.zi ./ 1000.0
     # update lateral subsurface flow domain (kinematic wave)
@@ -417,6 +467,8 @@ through BMI, to couple the SBM model to an external groundwater model.
 function update_until_recharge(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
     @unpack lateral, vertical, network, clock, config = model
 
+    do_water_demand = haskey(config.model, "water_demand")
+
     inds_riv = network.index_river
 
     # extract water levels h_av [m] from the land and river domains
@@ -438,6 +490,12 @@ function update_until_recharge(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:Sb
         )
     end
 
+    # optional water demand and allocation
+    if do_water_demand
+        update_water_demand(vertical)
+        update_water_allocation(model)
+    end
+
     # update vertical sbm concept until recharge [mm] to the saturated store
     update_until_recharge(vertical, config)
 
@@ -481,8 +539,7 @@ function update_total_water_storage(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,
     update_total_water_storage(
         vertical,
         network.index_river,
-        network.land.xl,
-        network.land.yl,
+        network.land.area,
         lateral.river,
         lateral.land,
     )
diff --git a/src/states.jl b/src/states.jl
index 9af01cc0e..d33427886 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -97,6 +97,10 @@ function extract_required_states(config::Config)
     do_lakes = get(config.model, "lakes", false)::Bool
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_floodplains = get(config.model, "floodplain_1d", false)::Bool
+    do_paddy = false
+    if haskey(config.model, "water_demand")
+        do_paddy = get(config.model.water_demand, "paddy", false)::Bool
+    end 
 
     # Extract required stated based on model configuration file
     vertical_states = get_vertical_states(model_type; snow = do_snow, glacier = do_glaciers)
@@ -165,6 +169,9 @@ function extract_required_states(config::Config)
     lake_states = do_lakes ? (:waterlevel,) : nothing
     reservoir_states = do_reservoirs ? (:volume,) : nothing
 
+    # Paddy states
+    paddy_states = do_paddy ? (:h,) : nothing
+
     # Build required states in a tuple, similar to the keys in the output of
     # `ncnames(config.state)`
     required_states = ()
@@ -198,7 +205,12 @@ function extract_required_states(config::Config)
         (:lateral, :river, :reservoir),
         reservoir_states,
     )
-
+    # Add paddy states to dict
+    required_states = add_to_required_states(
+        required_states,
+        (:vertical, :paddy),
+        paddy_states,
+    )
     return required_states
 end
 
diff --git a/src/utils.jl b/src/utils.jl
index 57563b7b6..1f557b1cb 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -810,7 +810,8 @@ end
 function initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile)
     for i in eachindex(ssf.ssf)
         ssf.kh[i] = kh_layered_profile(sbm, ssf.khfrac[i], i, ksat_profile)
-        ssf.ssf[i] = ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.slope[i] * ssf.dw[i]
+        ssf.ssf[i] =
+            ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.slope[i] * ssf.dw[i]
         kh_max = 0.0
         for j = 1:sbm.nlayers[i]
             if j <= sbm.nlayers_kv[i]
@@ -826,3 +827,14 @@ function initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile)
         ssf.ssfmax[i] = kh_max * ssf.slope[i]
     end
 end
+
+"""
+    bounded_divide(x, y; max = 1.0, default = 0.0)
+
+Return the division of `x` by `y`, bounded by a maximum value `max`, when `y` > 0.0.
+Otherwise return a `default` value.
+"""
+function bounded_divide(x, y; max = 1.0, default = 0.0)
+    z = y > 0.0 ? min(x / y, max) : default
+    return z
+end
diff --git a/src/vertical_process.jl b/src/vertical_process.jl
index 9801b1e2e..032a0fc6c 100644
--- a/src/vertical_process.jl
+++ b/src/vertical_process.jl
@@ -120,96 +120,78 @@ function rainfall_interception_modrut(
 end
 
 """
-    acttransp_unsat_sbm(rootingdepth, ustorelayerdepth, sumlayer, restpotevap, sum_actevapustore, c, usl, theta_s, theta_r, hb, ust::Bool = false)
+    vwc_brooks_corey(h, hb, theta_s, theta_r, c)
 
-Compute actual transpiration for unsaturated zone.
-If `ust` is `true`, the whole unsaturated store is available for transpiration.
-
-# Arguments
-- `rootingdepth`
-- `ustorelayerdepth`
-- `sumlayer` (depth (z) of upper boundary unsaturated layer)
-- `restpotevap` (remaining evaporation)
-- `sum_actevapustore` (cumulative actual transpiration (more than one unsaturated layers))
-- `c` (Brooks-Corey coefficient)
-- `usl` (thickness of unsaturated zone)
-- `theta_s`
-- `theta_r`
-- `hb` (air entry pressure)
-- `ust`
-
-# Output
-- `ustorelayerdepth`
-- `sum_actevapustore`
-- `restpotevap`
+Volumetric water content based on the Brooks-Corey soil hydraulic model.
 """
-function acttransp_unsat_sbm(
-    rootingdepth,
-    ustorelayerdepth,
-    sumlayer,
-    restpotevap,
-    sum_actevapustore,
-    c,
-    usl,
-    theta_s,
-    theta_r,
-    hb,
-    ust::Bool = false,
-)
-
-    # AvailCap is fraction of unsat zone containing roots
-    if ust
-        availcap = ustorelayerdepth * 0.99
+function vwc_brooks_corey(h, hb, theta_s, theta_r, c)
+    if h < hb
+        par_lambda = 2.0 / (c - 3.0)
+        vwc = (theta_s - theta_r) * pow(hb / h, par_lambda) + theta_r
     else
-        if usl > 0.0
-            availcap = min(1.0, max(0.0, (rootingdepth - sumlayer) / usl))
-        else
-            availcap = 0.0
-        end
+        vwc = theta_s
     end
+    return vwc
+end
 
-    maxextr = availcap * ustorelayerdepth
-
-    # Next step is to make use of the Feddes curve in order to decrease actevapustore when soil moisture values
-    # occur above or below ideal plant growing conditions (see also Feddes et al., 1978). h1-h4 values are
-    # actually negative, but all values are made positive for simplicity.
-    h1 = hb  # cm (air entry pressure)
-    h2 = 100.0  # cm (pF 2 for field capacity)
-    h3 = 400.0  # cm (pF 3, critical pF value)
-    h4 = 15849.0  # cm (pF 4.2, wilting point)
+"""
+    head_brooks_corey(vwc, theta_s, theta_r, c, hb)
 
-    # According to Brooks-Corey
+Soil water pressure head based on the Brooks-Corey soil hydraulic model.
+"""
+function head_brooks_corey(vwc, theta_s, theta_r, c, hb)
     par_lambda = 2.0 / (c - 3.0)
-    if usl > 0.0
-        vwc = ustorelayerdepth / usl
-    else
-        vwc = 0.0
-    end
-    vwc = max(vwc, 0.0000001)
-    head = hb / (pow(((vwc) / (theta_s - theta_r)), (1.0 / par_lambda)))  # Note that in the original formula, thetaR is extracted from vwc, but thetaR is not part of the numerical vwc calculation
-    head = max(head, hb)
-
-    # Transform h to a reduction coefficient value according to Feddes et al. (1978).
-    # For now: no reduction for head < h2 until following improvement (todo):
-    #       - reduction only applied to crops
-    if head <= h1
-        alpha = 1.0
-    elseif head >= h4
-        alpha = 0.0
-    elseif (head < h2) && (head > h1)
-        alpha = 1.0
-    elseif (head > h3) && (head < h4)
-        alpha = 1.0 - (head - h3) / (h4 - h3)
+    # Note that in the original formula, theta_r is extracted from vwc, but theta_r is not part of the numerical vwc calculation
+    h = hb / (pow(((vwc) / (theta_s - theta_r)), (1.0 / par_lambda)))
+    h = min(h, hb)
+    return h
+end
+
+"""
+    feddes_h3(h3_high, h3_low, tpot, Δt)
+
+Return soil water pressure head `h3` of Feddes root water uptake reduction function.
+"""
+function feddes_h3(h3_high, h3_low, tpot, Δt)
+    # value of h3 is a function of potential transpiration [mm/d]
+    tpot_daily = tpot * (basetimestep / Δt)
+    if (tpot_daily >= 0.0) && (tpot_daily <= 1.0)
+        h3 = h3_low
+    elseif (tpot_daily > 1.0) && (tpot_daily < 5.0)
+        h3 = h3_high + ((h3_low - h3_high) * (5.0 - tpot_daily)) / (5.0 - 1.0)
     else
-        alpha = 1.0
+        h3 = h3_high
     end
+    return h3
+end
 
-    actevapustore = (min(maxextr, restpotevap, ustorelayerdepth)) * alpha
-    ustorelayerdepth = ustorelayerdepth - actevapustore
-    restpotevap = restpotevap - actevapustore
-    sum_actevapustore = actevapustore + sum_actevapustore
+"""
+    rwu_reduction_feddes(h, h1, h2, h3, h4, alpha_h1)
 
-    return ustorelayerdepth, sum_actevapustore, restpotevap
+Root water uptake reduction factor based on Feddes.
+"""
+function rwu_reduction_feddes(h, h1, h2, h3, h4, alpha_h1)
+    # root water uptake reduction coefficient alpha (see also Feddes et al., 1978)
+    if alpha_h1 == 0.0
+        if (h <= h4) || (h > h1)
+            alpha = 0.0
+        elseif (h > h2) && (h <= h1)
+            alpha = (h - h1) / (h2 - h1)
+        elseif (h >= h3) && (h <= h2)
+            alpha = 1.0
+        elseif (h >= h4) && (h < h3)
+            alpha = (h - h4) / (h3 - h4)
+        end
+    else
+        if h <= h4
+            alpha = 0.0
+        elseif h >= h3
+            alpha = 1.0
+        elseif (h >= h4) && (h < h3)
+            alpha = (h - h4) / (h3 - h4)
+        end
+    end
+    return alpha
 end
 
 """
@@ -259,7 +241,13 @@ function infiltration(
 
     infiltexcess = (soilinf - max_infiltsoil) + (pathinf - max_infiltpath)
 
-    return infiltsoilpath, infiltsoil, infiltpath, soilinf, pathinf, infiltexcess
+    return infiltsoilpath,
+    infiltsoil,
+    infiltpath,
+    soilinf,
+    pathinf,
+    infiltexcess,
+    soilinfredu
 end
 
 """
diff --git a/src/water_demand.jl b/src/water_demand.jl
new file mode 100644
index 000000000..32e364e0c
--- /dev/null
+++ b/src/water_demand.jl
@@ -0,0 +1,626 @@
+@get_units @exchange @grid_type @grid_location @with_kw struct NonIrrigationDemand{T}
+    demand_gross::Vector{T}                 # gross water demand [mm Δt⁻¹]
+    demand_net::Vector{T}                   # net water demand [mm Δt⁻¹]
+    returnflow_fraction::Vector{T} | "-"    # return flow fraction [-]
+    returnflow::Vector{T}                   # return flow [mm Δt⁻¹]
+end
+
+@get_units @exchange @grid_type @grid_location @with_kw struct NonPaddy{T}
+    demand_gross::Vector{T}                     # irrigation gross demand [mm Δt⁻¹] 
+    irrigation_efficiency::Vector{T} | "-"      # irrigation efficiency [-]
+    maximum_irrigation_rate::Vector{T}          # maximum irrigation depth [mm Δt⁻¹]
+    irrigation_areas::Vector{Bool} | "-"        # irrigation areas [-]
+    irrigation_trigger::Vector{Bool} | "-"      # irrigation on or off [-]
+end
+
+@get_units @exchange @grid_type @grid_location @with_kw struct Paddy{T}
+    demand_gross::Vector{T}                     # irrigation gross demand [mm Δt⁻¹] 
+    irrigation_efficiency::Vector{T} | "-"      # irrigation efficiency [-]
+    maximum_irrigation_rate::Vector{T}          # maximum irrigation depth [mm Δt⁻¹]
+    irrigation_areas::Vector{Bool} | "-"        # irrigation areas [-]
+    irrigation_trigger::Vector{Bool} | "-"      # irrigation on or off [-]
+    h_min::Vector{T} | "mm"                     # minimum required water depth in the irrigated rice field [mm]
+    h_opt::Vector{T} | "mm"                     # optimal water depth in the irrigated rice fields [mm]
+    h_max::Vector{T} | "mm"                     # water depth when rice field starts spilling water (overflow) [mm]
+    h::Vector{T} | "mm"                         # actual water depth in rice field [mm]
+end
+
+@get_units @exchange @grid_type @grid_location @with_kw struct AllocationRiver{T}
+    act_surfacewater_abst::Vector{T}                    # actual surface water abstraction [mm Δt⁻¹]
+    act_surfacewater_abst_vol::Vector{T} | "m3 dt-1"    # actual surface water abstraction [m³ Δt⁻¹]
+    available_surfacewater::Vector{T} | "m3"            # available surface water [m³]
+    nonirri_returnflow::Vector{T}                       # return flow from non irrigation [mm Δt⁻¹] 
+end
+
+@get_units @exchange @grid_type @grid_location @with_kw struct AllocationLand{T}
+    irri_demand_gross::Vector{T}                        # irrigation gross demand [mm Δt⁻¹]
+    nonirri_demand_gross::Vector{T}                     # non-irrigation gross demand [mm Δt⁻¹]
+    total_gross_demand::Vector{T}                       # total gross demand [mm Δt⁻¹]
+    frac_sw_used::Vector{T} | "-"                       # fraction surface water used [-]
+    areas::Vector{Int} | "-"                            # allocation areas [-]
+    surfacewater_demand::Vector{T}                      # demand from surface water [mm Δt⁻¹]
+    surfacewater_alloc::Vector{T}                       # allocation from surface water [mm Δt⁻¹]
+    act_groundwater_abst::Vector{T}                     # actual groundwater abstraction [mm Δt⁻¹]
+    act_groundwater_abst_vol::Vector{T} | "m3 dt-1"     # actual groundwater abstraction [m³ Δt⁻¹]
+    available_groundwater::Vector{T} | "m3"             # available groundwater [m³]
+    groundwater_demand::Vector{T}                       # demand from groundwater [mm Δt⁻¹]
+    groundwater_alloc::Vector{T}                        # allocation from groundwater [mm Δt⁻¹]
+    irri_alloc::Vector{T}                               # allocated water for irrigation [mm Δt⁻¹]
+    nonirri_alloc::Vector{T}                            # allocated water for non-irrigation [mm Δt⁻¹]
+    total_alloc::Vector{T}                              # total allocated water [mm Δt⁻¹]
+    nonirri_returnflow::Vector{T}                       # return flow from non irrigation [mm Δt⁻¹]
+end
+
+"Return return flow fraction based on gross water demand `demand_gross` and net water demand `demand_net`"
+function returnflow_fraction(demand_gross, demand_net)
+    fraction = bounded_divide(demand_net, demand_gross)
+    returnflow_fraction = 1.0 - fraction
+    return returnflow_fraction
+end
+
+"Initialize water demand for the domestic sector"
+function initialize_domestic_demand(nc, config, inds, dt)
+    demand_gross =
+        ncread(
+            nc,
+            config,
+            "vertical.domestic.demand_gross";
+            sel = inds,
+            defaults = 0.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+    demand_net =
+        ncread(
+            nc,
+            config,
+            "vertical.domestic.demand_net";
+            sel = inds,
+            defaults = 0.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+    n = length(inds)
+    returnflow_f = returnflow_fraction.(demand_gross, demand_net)
+
+    domestic = NonIrrigationDemand{Float}(
+        demand_gross = demand_gross,
+        demand_net = demand_net,
+        returnflow_fraction = returnflow_f,
+        returnflow = fill(Float(0), n),
+    )
+
+    return domestic
+end
+
+"Initialize water demand for the industry sector"
+function initialize_industry_demand(nc, config, inds, dt)
+    demand_gross =
+        ncread(
+            nc,
+            config,
+            "vertical.industry.demand_gross";
+            sel = inds,
+            defaults = 0.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+    demand_net =
+        ncread(
+            nc,
+            config,
+            "vertical.industry.demand_net";
+            sel = inds,
+            defaults = 0.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+    n = length(inds)
+    returnflow_f = returnflow_fraction.(demand_gross, demand_net)
+
+    industry = NonIrrigationDemand{Float}(
+        demand_gross = demand_gross,
+        demand_net = demand_net,
+        returnflow_fraction = returnflow_f,
+        returnflow = fill(Float(0), n),
+    )
+
+    return industry
+end
+
+"Initialize water demand for the livestock sector"
+function initialize_livestock_demand(nc, config, inds, dt)
+    demand_gross =
+        ncread(
+            nc,
+            config,
+            "vertical.livestock.demand_gross";
+            sel = inds,
+            defaults = 0.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+    demand_net =
+        ncread(
+            nc,
+            config,
+            "vertical.livestock.demand_net";
+            sel = inds,
+            defaults = 0.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+    n = length(inds)
+    returnflow_f = returnflow_fraction.(demand_gross, demand_net)
+
+    livestock = NonIrrigationDemand{Float}(
+        demand_gross = demand_gross,
+        demand_net = demand_net,
+        returnflow_fraction = returnflow_f,
+        returnflow = fill(Float(0), n),
+    )
+    return livestock
+end
+
+"Initialize paddy (rice) fields for water demand and irrigation computations"
+function initialize_paddy(nc, config, inds, dt)
+    h_min = ncread(
+        nc,
+        config,
+        "vertical.paddy.h_min";
+        sel = inds,
+        defaults = 20.0,
+        type = Float,
+    )
+    h_opt = ncread(
+        nc,
+        config,
+        "vertical.paddy.h_opt";
+        sel = inds,
+        defaults = 50.0,
+        type = Float,
+    )
+    h_max = ncread(
+        nc,
+        config,
+        "vertical.paddy.h_max";
+        sel = inds,
+        defaults = 80.0,
+        type = Float,
+    )
+    efficiency = ncread(
+        nc,
+        config,
+        "vertical.paddy.irrigation_efficiency";
+        sel = inds,
+        defaults = 1.0,
+        type = Float,
+    )
+    areas = ncread(
+        nc,
+        config,
+        "vertical.paddy.irrigation_areas";
+        sel = inds,
+        optional = false,
+        type = Bool,
+    )
+    irrigation_trigger = ncread(
+        nc,
+        config,
+        "vertical.paddy.irrigation_trigger";
+        sel = inds,
+        optional = false,
+        type = Bool,
+    )
+    max_irri_rate =
+        ncread(
+            nc,
+            config,
+            "vertical.paddy.maximum_irrigation_rate";
+            sel = inds,
+            defaults = 25.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+
+    paddy = Paddy{Float}(
+        demand_gross = fill(mv, length(inds)),
+        irrigation_efficiency = efficiency,
+        maximum_irrigation_rate = max_irri_rate,
+        irrigation_trigger = irrigation_trigger,
+        h_min = h_min,
+        h_max = h_max,
+        h_opt = h_opt,
+        irrigation_areas = areas,
+        h = fill(0.0, length(inds)),
+    )
+    return paddy
+end
+
+"Initialize crop (non paddy) fields for water demand and irrigation computations"
+function initialize_nonpaddy(nc, config, inds, dt)
+    efficiency = ncread(
+        nc,
+        config,
+        "vertical.nonpaddy.irrigation_efficiency";
+        sel = inds,
+        defaults = 1.0,
+        type = Float,
+    )
+    areas = ncread(
+        nc,
+        config,
+        "vertical.nonpaddy.irrigation_areas";
+        sel = inds,
+        defaults = 1,
+        optional = false,
+        type = Int,
+    )
+    irrigation_trigger = ncread(
+        nc,
+        config,
+        "vertical.nonpaddy.irrigation_trigger";
+        sel = inds,
+        defaults = 1,
+        optional = false,
+        type = Bool,
+    )
+    max_irri_rate =
+        ncread(
+            nc,
+            config,
+            "vertical.nonpaddy.maximum_irrigation_rate";
+            sel = inds,
+            defaults = 25.0,
+            type = Float,
+        ) .* (dt / basetimestep)
+
+    nonpaddy = NonPaddy{Float}(
+        demand_gross = fill(mv, length(inds)),
+        maximum_irrigation_rate = max_irri_rate,
+        irrigation_efficiency = efficiency,
+        irrigation_areas = areas,
+        irrigation_trigger = irrigation_trigger,
+    )
+
+    return nonpaddy
+end
+
+"Initialize water allocation for the river domain"
+function initialize_allocation_river(n)
+    allocation = AllocationRiver(
+        act_surfacewater_abst = zeros(Float, n),
+        act_surfacewater_abst_vol = zeros(Float, n),
+        available_surfacewater = zeros(Float, n),
+        nonirri_returnflow = zeros(Float, n),
+    )
+    return allocation
+end
+
+"Initialize water allocation for the land domain (`vertical`)"
+function initialize_allocation_land(nc, config, inds)
+    frac_sw_used = ncread(
+        nc,
+        config,
+        "vertical.allocation.frac_sw_used";
+        sel = inds,
+        defaults = 1,
+        type = Float,
+    )
+    areas = ncread(
+        nc,
+        config,
+        "vertical.allocation.areas";
+        sel = inds,
+        defaults = 1,
+        type = Int,
+    )
+
+    n = length(inds)
+
+    allocation = AllocationLand(
+        irri_demand_gross = zeros(Float, n),
+        nonirri_demand_gross = zeros(Float, n),
+        total_gross_demand = zeros(Float, n),
+        frac_sw_used = frac_sw_used,
+        areas = areas,
+        surfacewater_demand = zeros(Float, n),
+        surfacewater_alloc = zeros(Float, n),
+        act_groundwater_abst = zeros(Float, n),
+        act_groundwater_abst_vol = zeros(Float, n),
+        available_groundwater = zeros(Float, n),
+        groundwater_demand = zeros(Float, n),
+        groundwater_alloc = zeros(Float, n),
+        irri_alloc = zeros(Float, n),
+        nonirri_alloc = zeros(Float, n),
+        total_alloc = zeros(Float, n),
+        nonirri_returnflow = zeros(Float, n),
+    )
+    return allocation
+end
+
+"Return non-irrigation gross demand and update returnflow fraction"
+function update_non_irrigation_demand(non_irri::NonIrrigationDemand, i)
+    non_irri.returnflow_fraction[i] =
+        returnflow_fraction(non_irri.demand_gross[i], non_irri.demand_net[i])
+    return non_irri.demand_gross[i]
+end
+
+# return zero (gross demand) if non-irrigation sector is not defined
+update_non_irrigation_demand(non_irri::Nothing, i) = 0.0
+
+"Update water allocation for river and land domains based on local surface water (river) availability."
+function surface_water_allocation_local(land, river, network)
+    # maps from the land domain to the internal river domain (linear index), excluding water bodies
+    index_river = network.land.index_river_wb
+    for i in eachindex(land.allocation.surfacewater_demand)
+        if index_river[i] > 0.0
+            # the available volume is limited by a fixed scaling factor of 0.8 to prevent
+            # rivers completely drying out. check for abstraction through inflow (external
+            # negative inflow) and adjust available volume.
+            if river.inflow[index_river[i]] < 0.0
+                inflow = river.inflow[index_river[i]] * land.dt
+                available_volume = max(river.volume[index_river[i]] * 0.80 + inflow, 0.0)
+            else
+                available_volume = river.volume[index_river[i]] * 0.80
+            end
+            # satisfy surface water demand with available local river volume
+            surfacewater_demand_vol =
+                land.allocation.surfacewater_demand[i] * 0.001 * network.land.area[i]
+            abstraction_vol = min(surfacewater_demand_vol, available_volume)
+            river.allocation.act_surfacewater_abst_vol[index_river[i]] = abstraction_vol
+            # remaining available surface water and demand 
+            river.allocation.available_surfacewater[index_river[i]] =
+                max(available_volume - abstraction_vol, 0.0)
+            abstraction = (abstraction_vol / network.land.area[i]) * 1000.0
+            land.allocation.surfacewater_demand[i] =
+                max(land.allocation.surfacewater_demand[i] - abstraction, 0.0)
+            # update actual abstraction from river and surface water allocation (land cell)
+            river.allocation.act_surfacewater_abst[index_river[i]] = abstraction
+            land.allocation.surfacewater_alloc[i] = abstraction
+        end
+    end
+end
+
+"Update water allocation for river and land domains based on surface water (river) availability for allocation areas."
+function surface_water_allocation_area(land, river, network)
+    inds_river = network.river.indices_allocation_areas
+    inds_land = network.land.indices_allocation_areas
+    res_index = network.river.reservoir_index
+    lake_index = network.river.lake_index
+
+    m = length(inds_river)
+    # loop over allocation areas
+    for i = 1:m
+        # surface water demand (allocation area)
+        sw_demand_vol = 0.0
+        for j in inds_land[i]
+            sw_demand_vol +=
+                land.allocation.surfacewater_demand[j] * 0.001 * network.land.area[j]
+        end
+        # surface water availability (allocation area)
+        sw_available = 0.0
+        for j in inds_river[i]
+            # for reservoir locations use reservoir volume
+            if res_index[j] > 0
+                k = res_index[j]
+                river.allocation.available_surfacewater[j] =
+                    river.reservoir.volume[k] * 0.98 # limit available reservoir volume
+                sw_available += river.allocation.available_surfacewater[j]
+                # for lake locations use lake volume
+            elseif lake_index[j] > 0
+                k = lake_index[j]
+                river.allocation.available_surfacewater[j] = river.lake.storage[k] * 0.98 # limit available lake volume
+                sw_available += river.allocation.available_surfacewater[j]
+                # river volume
+            else
+                sw_available += river.allocation.available_surfacewater[j]
+            end
+        end
+        # total actual surface water abstraction [m3] in an allocation area, minimum of
+        # available surface water and demand in an allocation area.
+        sw_abstraction = min(sw_available, sw_demand_vol)
+
+        # fraction of available surface water that can be abstracted at allocation area
+        # level
+        frac_abstract_sw = bounded_divide(sw_abstraction, sw_available)
+        # fraction of water demand that can be satisfied by available surface water at
+        # allocation area level. 
+        frac_allocate_sw = bounded_divide(sw_abstraction, sw_demand_vol)
+
+        # water abstracted from surface water at each river cell (including reservoir and
+        # lake locations).
+        for j in inds_river[i]
+            river.allocation.act_surfacewater_abst_vol[j] +=
+                frac_abstract_sw * river.allocation.available_surfacewater[j]
+            river.allocation.act_surfacewater_abst[j] =
+                (river.allocation.act_surfacewater_abst_vol[j] / network.river.area[j]) *
+                1000.0
+        end
+
+        # water allocated to each land cell.
+        for j in inds_land[i]
+            land.allocation.surfacewater_alloc[j] +=
+                frac_allocate_sw * land.allocation.surfacewater_demand[j]
+        end
+    end
+end
+
+"Update water allocation for subsurface domain based on local groundwater availability."
+function groundwater_allocation_local(land, groundwater_volume, network)
+    for i in eachindex(land.allocation.groundwater_demand)
+        # groundwater demand based on allocation from surface water.
+        land.allocation.groundwater_demand[i] = max(
+            land.allocation.total_gross_demand[i] - land.allocation.surfacewater_alloc[i],
+            0.0,
+        )
+        # land index excluding water bodies
+        if network.index_wb[i]
+            # satisfy groundwater demand with available local groundwater volume
+            groundwater_demand_vol =
+                land.allocation.groundwater_demand[i] * 0.001 * network.area[i]
+            available_volume = groundwater_volume[i] * 0.75 # limit available groundwater volume
+            abstraction_vol = min(groundwater_demand_vol, available_volume)
+            land.allocation.act_groundwater_abst_vol[i] = abstraction_vol
+            # remaining available groundwater and demand 
+            land.allocation.available_groundwater[i] =
+                max(available_volume - abstraction_vol, 0.0)
+            abstraction = (abstraction_vol / network.area[i]) * 1000.0
+            land.allocation.groundwater_demand[i] =
+                max(land.allocation.groundwater_demand[i] - abstraction, 0.0)
+            # update actual abstraction from groundwater and groundwater allocation (land cell)
+            land.allocation.act_groundwater_abst[i] = abstraction
+            land.allocation.groundwater_alloc[i] = abstraction
+        end
+    end
+end
+
+"Update water allocation for subsurface domain based on groundwater availability for allocation areas."
+function groundwater_allocation_area(land, network)
+    inds_river = network.river.indices_allocation_areas
+    inds_land = network.land.indices_allocation_areas
+    m = length(inds_river)
+    # loop over allocation areas
+    for i = 1:m
+        # groundwater demand and availability (allocation area)
+        gw_demand_vol = 0.0
+        gw_available = 0.0
+        for j in inds_land[i]
+            gw_demand_vol +=
+                land.allocation.groundwater_demand[j] * 0.001 * network.land.area[j]
+            gw_available += land.allocation.available_groundwater[j]
+        end
+        # total actual groundwater abstraction [m3] in an allocation area, minimum of
+        # available  groundwater and demand in an allocation area.
+        gw_abstraction = min(gw_available, gw_demand_vol)
+
+        # fraction of available groundwater that can be abstracted at allocation area level
+        frac_abstract_gw = bounded_divide(gw_abstraction, gw_available)
+        # fraction of water demand that can be satisfied by available groundwater at
+        # allocation area level.
+        frac_allocate_gw = bounded_divide(gw_abstraction, gw_demand_vol)
+
+        # water abstracted from groundwater and allocated.
+        for j in inds_land[i]
+            land.allocation.act_groundwater_abst_vol[j] +=
+                frac_abstract_gw * land.allocation.available_groundwater[j]
+            land.allocation.act_groundwater_abst[j] =
+                1000.0 *
+                (land.allocation.act_groundwater_abst_vol[j] / network.land.area[j])
+            land.allocation.groundwater_alloc[j] +=
+                frac_allocate_gw * land.allocation.groundwater_demand[j]
+        end
+    end
+end
+
+"Return and update non-irrigation sector (domestic, livestock, industry) return flow"
+function return_flow(non_irri::NonIrrigationDemand, allocation)
+    for i in eachindex(non_irri.returnflow)
+        frac = bounded_divide(non_irri.demand_gross[i], allocation.nonirri_demand_gross[i])
+        allocate = frac * allocation.nonirri_alloc[i]
+        non_irri.returnflow[i] = non_irri.returnflow_fraction[i] * allocate
+    end
+    return non_irri.returnflow
+end
+
+# return zero (return flow) if non-irrigation sector is not defined
+return_flow(non_irri::Nothing, allocation) = 0.0
+
+groundwater_volume(gw::LateralSSF) = gw.volume
+groundwater_volume(gw) = gw.flow.aquifer.volume
+
+"""
+    update_water_allocation(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:Union{SbmModel, SbmGwfModel}}
+
+Update water allocation for model type `sbm` or `sbm_gwf` for a single timestep. First,
+surface water abstraction is computed to satisfy local water demand (non-irrigation and
+irrigation), and then updated (including lakes and reservoirs) to satisfy the remaining
+water demand for allocation areas. Then groundwater abstraction is computed to satisfy the
+remaining local water demand, and then updated to satisfy the remaining water demand for
+allocation areas. Finally, non-irrigation return flows are updated.
+"""
+function update_water_allocation(
+    model::Model{N,L,V,R,W,T},
+) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
+
+    @unpack network, lateral, vertical = model
+
+    river = lateral.river
+    index_river = network.land.index_river_wb
+    res_index_f = network.river.reservoir_index_f
+    lake_index_f = network.river.lake_index_f
+
+    vertical.allocation.surfacewater_alloc .= 0.0
+    river.allocation.act_surfacewater_abst .= 0.0
+    river.allocation.act_surfacewater_abst_vol .= 0.0
+    # total surface water demand for each land cell
+    @. vertical.allocation.surfacewater_demand =
+        vertical.allocation.frac_sw_used * vertical.allocation.nonirri_demand_gross +
+        vertical.allocation.frac_sw_used * vertical.allocation.irri_demand_gross
+
+    # local surface water demand and allocation (river, excluding reservoirs and lakes)
+    surface_water_allocation_local(vertical, river, network)
+    # surface water demand and allocation for areas
+    surface_water_allocation_area(vertical, river, network)
+
+    @. river.abstraction = river.allocation.act_surfacewater_abst_vol / vertical.dt
+
+    # for reservoir and lake locations set river abstraction at zero and abstract volume
+    # from reservoir and lake, including an update of lake waterlevel
+    if !isnothing(river.reservoir)
+        @. river.abstraction[res_index_f] = 0.0
+        @. river.reservoir.volume -= river.allocation.act_surfacewater_abst_vol[res_index_f]
+    elseif !isnothing(river.lake)
+        @. river.abstraction[lake_index_f] = 0.0
+        lakes = river.lake
+        @. lakes.storage -= river.allocation.act_surfacewater_abst_vol[lake_index_f]
+        @. lakes.waterlevel =
+            waterlevel(lakes.storfunc, lakes.area, lakes.storage, lakes.sh)
+    end
+
+    vertical.allocation.groundwater_alloc .= 0.0
+    vertical.allocation.act_groundwater_abst_vol .= 0.0
+    vertical.allocation.act_groundwater_abst .= 0.0
+    # local groundwater demand and allocation
+    groundwater_allocation_local(
+        vertical,
+        groundwater_volume(lateral.subsurface),
+        network.land,
+    )
+    # groundwater demand and allocation for areas
+    groundwater_allocation_area(vertical, network)
+
+    # irrigation allocation
+    for i in eachindex(vertical.allocation.total_alloc)
+        vertical.allocation.total_alloc[i] =
+            vertical.allocation.groundwater_alloc[i] +
+            vertical.allocation.surfacewater_alloc[i]
+        frac_irri = bounded_divide(
+            vertical.allocation.irri_demand_gross[i],
+            vertical.allocation.total_gross_demand[i],
+        )
+        vertical.allocation.irri_alloc[i] = frac_irri * vertical.allocation.total_alloc[i]
+        vertical.allocation.nonirri_alloc[i] =
+            vertical.allocation.total_alloc[i] - vertical.allocation.irri_alloc[i]
+    end
+
+    # non-irrigation return flows
+    returnflow_livestock = return_flow(vertical.livestock, vertical.allocation)
+    returnflow_domestic = return_flow(vertical.domestic, vertical.allocation)
+    returnflow_industry = return_flow(vertical.industry, vertical.allocation)
+
+    # map non-irrigation return flow to land and river water allocation
+    if (
+        !isnothing(vertical.livestock) ||
+        !isnothing(vertical.domestic) ||
+        !isnothing(vertical.industry)
+    )
+        @. vertical.allocation.nonirri_returnflow =
+            returnflow_livestock + returnflow_domestic + returnflow_industry
+
+        for i in eachindex(vertical.allocation.nonirri_returnflow)
+            if index_river[i] > 0.0
+                k = index_river[i]
+                river.allocation.nonirri_returnflow[k] =
+                    vertical.allocation.nonirri_returnflow[i]
+                vertical.allocation.nonirri_returnflow[i] = 0.0
+            else
+                vertical.allocation.nonirri_returnflow[i] =
+                    vertical.allocation.nonirri_returnflow[i]
+            end
+        end
+    end
+end


Test patch:
diff --git a/build/create_binaries/download_test_data.jl b/build/create_binaries/download_test_data.jl
index d61310f3b..0ed330ed6 100644
--- a/build/create_binaries/download_test_data.jl
+++ b/build/create_binaries/download_test_data.jl
@@ -51,3 +51,9 @@ instates_sbm_gw_path =
 lake_sh_1_path = testdata(v"0.2.1", "lake_sh_1.csv", "lake_sh_1.csv")
 lake_sh_2_path = testdata(v"0.2.1", "lake_sh_2.csv", "lake_sh_2.csv")
 lake_hq_2_path = testdata(v"0.2.1", "lake_hq_2.csv", "lake_hq_2.csv")
+forcing_calendar_noleap_path =
+    testdata(v"0.2.8", "forcing-calendar-noleap.nc", "forcing-calendar-noleap.nc")
+forcing_piave_path = testdata(v"0.2.9", "inmaps-era5-2010-piave.nc", "forcing-piave.nc")
+staticmaps_piave_path = testdata(v"0.2.9", "staticmaps-piave.nc", "staticmaps-piave.nc")
+instates_piave_path = testdata(v"0.2.9", "instates-piave.nc", "instates-piave.nc")
+instates_piave_gwf_path = testdata(v"0.2.9", "instates-piave-gwf.nc", "instates-piave-gwf.nc")
\ No newline at end of file
diff --git a/server/test/client.jl b/server/test/client.jl
index d5eaa70a9..36e6b5dc5 100644
--- a/server/test/client.jl
+++ b/server/test/client.jl
@@ -32,8 +32,8 @@ end
 
 @testset "model information functions" begin
     @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm")
-    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 193)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 193)
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 207)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 207)
     to_check = [
         "vertical.nlayers",
         "vertical.theta_r",
@@ -66,7 +66,7 @@ vwc_1_size = 0
     vwc_1_size = Int(vwc_1_nbytes / vwc_1_itemsize)
     @test request((fn = "get_var_grid", name = "lateral.river.h")) == Dict("var_grid" => 3)
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
-    @test mean(request(msg)["value"]) ≈ 277.8362161218515
+    @test mean(request(msg)["value"]) ≈ 277.3620724821974
     msg = (fn = "get_value_ptr", name = "vertical.theta_s")
     @test mean(request(msg)["value_ptr"]) ≈ 0.4409211971535584
     msg = (
@@ -76,7 +76,7 @@ vwc_1_size = 0
         inds = [1, 5, 10],
     )
     @test request(msg)["value_at_indices"] ≈
-          [2.196562540141252f0, 2.6853997213000866f0, 3.4814142040324985f0]
+          [2.198747900215207f0, 2.6880427720508515f0, 3.4848783702629564f0]
     msg = (fn = "set_value", name = "vertical.zi", src = fill(300.0, zi_size))
     @test request(msg) == Dict("status" => "OK")
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
diff --git a/test/bmi.jl b/test/bmi.jl
index db8c95462..ba0f3751b 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -20,8 +20,8 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "model information functions" begin
             @test BMI.get_component_name(model) == "sbm"
-            @test BMI.get_input_item_count(model) == 193
-            @test BMI.get_output_item_count(model) == 193
+            @test BMI.get_input_item_count(model) == 207
+            @test BMI.get_output_item_count(model) == 207
             to_check = [
                 "vertical.nlayers",
                 "vertical.theta_r",
@@ -57,7 +57,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test_throws ErrorException BMI.get_value_ptr(model, "vertical.")
             dest = zeros(Float, size(model.vertical.zi))
             BMI.get_value(model, "vertical.zi", dest)
-            @test mean(dest) ≈ 276.16325589542333
+            @test mean(dest) ≈ 276.1625022866973
             @test BMI.get_value_at_indices(
                 model,
                 "vertical.vwc[1]",
@@ -75,7 +75,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
                 "lateral.river.q",
                 zeros(Float, 3),
                 [1, 100, 5617],
-            ) ≈ [0.623325399343309, 5.227139951657074, 0.02794287432778194]
+            ) ≈ [0.623325399343309, 5.227139951657074, 0.027942874327781947]
             BMI.set_value(model, "vertical.zi", fill(300.0, length(model.vertical.zi)))
             @test mean(
                 BMI.get_value(model, "vertical.zi", zeros(Float, size(model.vertical.zi))),
@@ -156,8 +156,8 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             sbm = model.vertical
             @test sbm.interception[1] ≈ 0.32734913737568716f0
             @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
-            @test sbm.snow[1] ≈ 3.484789961176288f0
-            @test sbm.recharge[5] ≈ -0.0f0
+            @test sbm.snow[1] ≈ 3.4847899611762876f0
+            @test sbm.recharge[5] ≈ 0.0f0
             @test sbm.zi[5] ≈ 300.0f0
         end
 
@@ -180,7 +180,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             sub = model.lateral.subsurface
             @test sbm.interception[1] ≈ 0.32734913737568716f0
             @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
-            @test sbm.snow[1] ≈ 3.484789961176288f0
+            @test sbm.snow[1] ≈ 3.4847899611762876f0
             @test sbm.recharge[5] ≈ 0.0f0
             @test sbm.zi[5] ≈ 250.0f0
             @test sub.zi[5] ≈ 0.25f0
diff --git a/test/groundwater.jl b/test/groundwater.jl
index 13d529ed1..5a59a3e20 100644
--- a/test/groundwater.jl
+++ b/test/groundwater.jl
@@ -46,6 +46,7 @@ function homogenous_aquifer(nrow, ncol)
         fill(0.1, ncell), # specific storage
         fill(1.0, ncell),  # storativity
         fill(0.0, connectivity.nconnection),  # conductance
+        fill(0.0, ncell), # total volume that can be released
     )
     unconf_aqf = Wflow.UnconfinedAquifer(
         [0.0, 7.5, 20.0],  # head
@@ -55,6 +56,7 @@ function homogenous_aquifer(nrow, ncol)
         fill(100.0, ncell),  # area
         fill(0.15, ncell),  # specific yield
         fill(0.0, connectivity.nconnection),  # conductance
+        fill(0.0, ncell), # total volume that can be released
         fill(3.0, ncell), # conductance reduction factor
     )
     return (connectivity, conf_aqf, unconf_aqf)
@@ -411,6 +413,7 @@ end
             fill(cellsize * cellsize, ncell),
             fill(specific_yield, ncell),
             fill(0.0, connectivity.nconnection),
+            fill(0.0, ncell),
             fill(gwf_f, ncell),
         )
         # constant head on left boundary, 0 at 0
@@ -435,7 +438,14 @@ end
         end
 
         head_analytical = [
-            transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta) for x in xc
+            transient_aquifer_1d(
+                x,
+                time,
+                conductivity,
+                specific_yield,
+                aquifer_length,
+                beta,
+            ) for x in xc
         ]
         difference = gwf.aquifer.head .- head_analytical
         # @test all(difference .< ?)  #TODO
@@ -471,6 +481,7 @@ end
             fill(cellsize * cellsize, ncell),
             fill(specific_yield, ncell),
             fill(0.0, connectivity.nconnection),
+            fill(0.0, ncell),
             fill(gwf_f, ncell),
         )
         # constant head on left boundary, 0 at 0
@@ -495,7 +506,14 @@ end
         end
 
         head_analytical = [
-            transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta) for x in xc
+            transient_aquifer_1d(
+                x,
+                time,
+                conductivity,
+                specific_yield,
+                aquifer_length,
+                beta,
+            ) for x in xc
         ]
         difference = gwf.aquifer.head .- head_analytical
         # @test all(difference .< ?)  #TODO
@@ -535,6 +553,7 @@ end
             fill(specific_storage, ncell),
             fill(storativity, ncell),
             fill(0.0, connectivity.nconnection), # conductance, to be set
+            fill(0.0, ncell), # total volume that can be released, to be set
         )
 
         cell_index = reshape(collect(range(1, ncell, step = 1)), shape)
diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl
index 3e9eb41ed..eb3c96763 100644
--- a/test/horizontal_process.jl
+++ b/test/horizontal_process.jl
@@ -210,6 +210,7 @@ end
         volume = fill(0.0, n),
         error = zeros(n),
         inflow = zeros(n),
+        abstraction = zeros(n),
         inflow_wb = zeros(n),
         inwater = zeros(n),
         dl = dl,
@@ -224,6 +225,7 @@ end
         reservoir = nothing,
         lake = nothing,
         floodplain = nothing,
+        allocation = nothing,
     )
 
     # run until steady state is reached
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index bc2446f92..aec71a06e 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -53,6 +53,7 @@ end
     )
 
     Wflow.update(lake, 1, 2500.0, 181, 86400.0)
+    @test Wflow.waterlevel(lake.storfunc, lake.area, lake.storage, lake.sh)[1] ≈ 19.672653848925634
     @test lake.outflow[1] ≈ 85.14292808113598
     @test lake.totaloutflow[1] ≈ 7.356348986210149e6
     @test lake.storage[1] ≈ 3.55111879238499e9
@@ -152,6 +153,7 @@ end
     )
 
     Wflow.update(lake, 1, 1500.0, 15, 86400.0)
+    @test Wflow.waterlevel(lake.storfunc, lake.area, lake.storage, lake.sh) ≈ [398.0] atol = 1e-2
     @test lake.outflow ≈ [1303.67476852] atol = 1e-2
     @test lake.totaloutflow ≈ [11.26375000e7] atol = 1e3
     @test lake.storage ≈ [4.293225e8] atol = 1e4
diff --git a/test/run.jl b/test/run.jl
index 084b11b93..45cd79750 100644
--- a/test/run.jl
+++ b/test/run.jl
@@ -3,9 +3,6 @@ using Wflow
 using Dates
 using NCDatasets
 
-tomlpath = joinpath(@__DIR__, "sbm_simple.toml")
-Wflow.run(tomlpath; silent = true)
-
 # test whether restarted runs get the same results as continuous ones, i.e. state is captured
 tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 config = Wflow.Config(tomlpath)
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 62fdfba99..d50ccb752 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -34,7 +34,7 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     @test row.Q_6136151 ≈ 0.007946301730645885f0
     @test row.Q_6136160 ≈ 3.927719795530719f0
     @test row.Q_6136202 ≈ 1.4162246003743886f0
-    @test row.recharge_1 ≈ -0.002257181032501202f0
+    @test row.recharge_1 ≈ -0.0020800523945940217f0
 end
 
 @testset "NetCDF scalar output" begin
@@ -79,7 +79,7 @@ end
     @test sbm.theta_r[50063] ≈ 0.15943120419979095f0
     @test mean(sbm.runoff) ≈ 0.04177459898728149f0
     @test mean(sbm.soilevap) ≈ 0.02122698830889417f0
-    @test mean(sbm.actevap) ≈ 0.33545623834952154f0
+    @test mean(sbm.actevap) ≈ 0.3353001180202587f0
     @test mean(sbm.actinfilt) ≈ 1.6444774688444848f0
     @test sbm.snow[5] ≈ 3.768513390588815f0
     @test mean(sbm.snow) ≈ 0.038019723676094325f0
@@ -94,10 +94,10 @@ model = Wflow.run_timestep(model)
     sbm = model.vertical
     @test sbm.theta_s[50063] ≈ 0.48755401372909546f0
     @test sbm.theta_r[50063] ≈ 0.15943120419979095f0
-    @test mean(sbm.net_runoff) ≈ 0.23710522236947756f0
-    @test mean(sbm.runoff) ≈ 0.23747368431143517f0
+    @test mean(sbm.net_runoff) ≈ 0.23734052031823816f0
+    @test mean(sbm.runoff) ≈ 0.23770898226019577f0
     @test mean(sbm.soilevap) ≈ 0.018750808322054897f0
-    @test mean(sbm.actevap) ≈ 0.27379330211866343f0
+    @test mean(sbm.actevap) ≈ 0.14545276216428166f0
     @test mean(sbm.actinfilt) ≈ 0.08863102527394363f0
     @test sbm.snow[5] ≈ 3.843412524052313f0
     @test mean(sbm.snow) ≈ 0.03461317061870949f0
@@ -107,15 +107,15 @@ end
 
 @testset "subsurface flow" begin
     ssf = model.lateral.subsurface.ssf
-    @test sum(ssf) ≈ 6.370399148012509f7
-    @test ssf[network.land.order[1]] ≈ 7.169036749244327f2
-    @test ssf[network.land.order[end-100]] ≈ 2333.801056570759f0
-    @test ssf[network.land.order[end]] ≈ 288.19428729403944f0
+    @test sum(ssf) ≈ 6.3761585406186976f7
+    @test ssf[network.land.order[1]] ≈ 718.2802566393531f0
+    @test ssf[network.land.order[end-100]] ≈ 2337.771227118579f0
+    @test ssf[network.land.order[end]] ≈ 288.19428729403984f0
 end
 
 @testset "overland flow" begin
     q = model.lateral.land.q_av
-    @test sum(q) ≈ 291.27639107427285f0
+    @test sum(q) ≈ 291.4923871784623f0
     @test q[26625] ≈ 0.0
     @test q[39308] ≈ 0.0
     @test q[network.land.order[end]] ≈ 1.0f-30
@@ -123,9 +123,9 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3622.7369292570543f0
-    @test q[1622] ≈ 0.0006497468064774366f0
-    @test q[43] ≈ 12.05767242907667f0
+    @test sum(q) ≈ 3625.0013368279815f0
+    @test q[1622] ≈ 0.0006503254947860838f0
+    @test q[43] ≈ 12.06416878694095f0
     @test q[network.river.order[end]] ≈ 0.039200124520463835f0
 end
 
@@ -166,10 +166,10 @@ end
 
 @testset "river flow at basin outlets and downstream of one pit" begin
     q = model.lateral.river.q_av
-    @test q[4009] ≈ 8.51907041734622f0 # pit/ outlet, CartesianIndex(141, 228)
+    @test q[4009] ≈ 8.543145028037452f0 # pit/ outlet, CartesianIndex(141, 228)
     @test q[4020] ≈ 0.006779014715290862f0 # downstream of pit 4009, CartesianIndex(141, 229)
-    @test q[2508] ≈ 150.28398167251638f0 # pit/ outlet
-    @test q[5808] ≈ 0.12419895007970105f0 # pit/ outlet
+    @test q[2508] ≈ 150.5640617045796f0 # pit/ outlet
+    @test q[5808] ≈ 0.12438899196040518f0 # pit/ outlet
 end
 
 # test changing forcing and cyclic LAI parameter
@@ -212,7 +212,7 @@ model = Wflow.run_timestep(model)
 
 @testset "river inflow (cyclic)" begin
     @test model.lateral.river.inflow[44] ≈ 0.75
-    @test model.lateral.river.q_av[44] ≈ 10.71846407068599
+    @test model.lateral.river.q_av[44] ≈ 10.723729440690567f0
 end
 
 # test fixed forcing (precipitation = 2.5)
@@ -248,14 +248,14 @@ model = Wflow.run_timestep(model)
 
 @testset "river flow and depth (local inertial)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3919.6025219496014f0
-    @test q[1622] ≈ 7.31010246736994f-5
-    @test q[43] ≈ 11.92153120707289f0
-    @test q[501] ≈ 3.5736389982451895f0
+    @test sum(q) ≈ 3922.0644366362544f0
+    @test q[1622] ≈ 7.315676375562105f-5
+    @test q[43] ≈ 11.92787156357907f0
+    @test q[501] ≈ 3.57855182713785f0
     h = model.lateral.river.h_av
     @test h[1622] ≈ 0.001987887644883981f0
     @test h[43] ≈ 0.4366415244811759f0
-    @test h[501] ≈ 0.057265962518284294f0
+    @test h[501] ≈ 0.057317706869865745f0
     q_channel = model.lateral.river.q_channel_av
     @test q ≈ q_channel
 end
@@ -272,20 +272,20 @@ model = Wflow.run_timestep(model)
 @testset "river and overland flow and depth (local inertial)" begin
     q = model.lateral.river.q_av
     @test sum(q) ≈ 2380.64389229669f0
-    @test q[1622] ≈ 7.322956970529551f-5
-    @test q[43] ≈ 5.361283165612762f0
-    @test q[501] ≈ 1.6021771576366957f0
+    @test q[1622] ≈ 7.328535246760549f-5
+    @test q[43] ≈ 5.3566292152594155f0
+    @test q[501] ≈ 1.6042388573126602f0
     h = model.lateral.river.h_av
     @test h[1622] ≈ 0.0019891342000364796f0
-    @test h[43] ≈ 0.3003008810153667f0
-    @test h[501] ≈ 0.031925992442532f0
+    @test h[43] ≈ 0.30026439683630496f0
+    @test h[501] ≈ 0.03195324587192846f0
     qx = model.lateral.land.qx
     qy = model.lateral.land.qy
-    @test qx[[26, 35, 631]] ≈ [0.18614776104106373f0, 0.029502872625766417f0, 0.0f0]
-    @test qy[[26, 35, 631]] ≈ [0.12757214437549858f0, 1.7212079599401755f0, 0.0f0]
+    @test qx[[26, 35, 631]] ≈ [0.1939736998417174f0, 0.026579954465883678f0, 0.0f0]
+    @test qy[[26, 35, 631]] ≈ [0.12906530420401777f0, 1.7225115950614904f0, 0.0f0]
     h = model.lateral.land.h
     @test h[[26, 35, 631]] ≈
-          [0.07361854999908582f0, 0.009155393111676267f0, 0.0007258741013439351f0]
+          [0.07367301172613304f0, 0.009139882310161706f0, 0.0007482998926237368f0]
 end
 Wflow.close_files(model, delete_output = false)
 
@@ -401,16 +401,16 @@ model = Wflow.run_timestep(model)
 
 @testset "river flow (local inertial) with floodplain schematization simulation" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3908.039208613999f0
-    @test q[1622] ≈ 7.310102468091527f-5
-    @test q[43] ≈ 11.921531207072922f0
-    @test q[501] ≈ 3.5061516913374717f0
-    @test q[5808] ≈ 0.0022258226497210145f0
+    @test sum(q) ≈ 3910.4728717811836f0
+    @test q[1622] ≈ 7.315676384849305f-5
+    @test q[43] ≈ 11.92787156357908f0
+    @test q[501] ≈ 3.510668846752431f0
+    @test q[5808] ≈ 0.002223993845806248f0
     h = model.lateral.river.h_av
     @test h[1622] ≈ 0.001987887580593841f0
     @test h[43] ≈ 0.436641524481545f0
-    @test h[501] ≈ 0.05665942153713204f0
-	@test h[5808] ≈ 0.005933025055544241f0
+    @test h[501] ≈ 0.05670770509802258f0
+    @test h[5808] ≈ 0.005929945681367346f0
 end
 
 # set boundary condition local inertial routing from netCDF file
@@ -422,15 +422,15 @@ model = Wflow.run_timestep(model)
 
 @testset "change boundary condition for local inertial routing (including floodplain)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3908.2493947503826f0
-    @test q[1622] ≈ 7.310183576829848f-5
-    @test q[43] ≈ 11.921531207085014f0
-    @test q[501] ≈ 3.5061513868779763f0
+    @test sum(q) ≈ 3910.683449719468f0
+    @test q[1622] ≈ 7.315757521099307f-5
+    @test q[43] ≈ 11.927871563591228f0
+    @test q[501] ≈ 3.5106678593721496f0
     @test q[5808] ≈ 0.060518234525259465f0
     h = model.lateral.river.h_av
     @test h[1622] ≈ 0.0019878952928530183f0
     @test h[43] ≈ 0.4366415249636809f0
-    @test h[501] ≈ 0.05665929962422606f0
+    @test h[501] ≈ 0.056707564314724804f0
     @test h[5808] ≈ 2.0000006940603936f0
 end
 Wflow.close_files(model, delete_output = false)
@@ -533,9 +533,9 @@ Wflow.close_files(model, delete_output = false)
     model = Wflow.run_timestep(model)
     @testset "river flow layered exponential profile" begin
         q = model.lateral.river.q_av
-        @test sum(q) ≈ 3158.1868156296678f0
+        @test sum(q) ≈ 3159.38300016008f0
         @test q[1622] ≈ 0.0005972577112819149f0
-        @test q[43] ≈ 10.013363662625276f0
+        @test q[43] ≈ 10.017642376280731f0
     end
 
     Wflow.close_files(model, delete_output = false)
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index c2d422f03..b984fed2b 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -14,7 +14,7 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
 
     @test row.time == DateTime("2000-06-01T00:00:00")
     @test row.Q_av ≈ 0.01620324716944374f0
-    @test row.head ≈ 1.8416896245327632f0
+    @test row.head ≈ 1.6471323360175287f0
 end
 
 @testset "first timestep" begin
@@ -24,7 +24,7 @@ end
     @test sbm.theta_s[1] ≈ 0.44999998807907104f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
-    @test sbm.transpiration[1] ≈ 0.6117526566330049f0
+    @test sbm.transpiration[1] ≈ 0.30587632831650247f0
 end
 
 # run the second timestep
@@ -35,36 +35,36 @@ model = Wflow.run_timestep(model)
     @test sbm.theta_s[1] ≈ 0.44999998807907104f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
-    @test sbm.transpiration[1] ≈ 1.0122634204681036f0
+    @test sbm.transpiration[4] ≈ 0.7000003898938235f0
 end
 
 @testset "overland flow (kinematic wave)" begin
     q = model.lateral.land.q_av
-    @test sum(q) ≈ 2.2298616f-7
+    @test sum(q) ≈ 2.229860508650628f-7
 end
 
 @testset "river domain (kinematic wave)" begin
     q = model.lateral.river.q_av
     river = model.lateral.river
-    @test sum(q) ≈ 0.034904078613089404f0
-    @test q[6] ≈ 0.007881734464420193f0
-    @test river.volume[6] ≈ 4.479866566743959f0
-    @test river.inwater[6] ≈ 0.00036640965826798314f0
-    @test q[13] ≈ 0.0005958036525897934f0
-    @test q[network.river.order[end]] ≈ 0.008405059891856978f0
+    @test sum(q) ≈ 0.035443370536496675f0
+    @test q[6] ≈ 0.008031554512314907f0
+    @test river.volume[6] ≈ 4.532124903256408f0
+    @test river.inwater[6] ≈ 0.0004073892212290558f0
+    @test q[13] ≈ 0.0006017024138583771f0
+    @test q[network.river.order[end]] ≈ 0.008559590281509943f0
 end
 
 @testset "groundwater" begin
     gw = model.lateral.subsurface
     @test gw.river.stage[1] ≈ 1.2123636929067039f0
     @test gw.flow.aquifer.head[17:21] ≈ [
-        1.288882128227083f0,
-        1.344816977641676f0,
+        1.2866380350225155f0,
+        1.3477853512604643f0,
         1.7999999523162842f0,
-        1.803682542252168f0,
-        1.4049539807162532f0,
+        1.6225103807809076f0,
+        1.4053590307668113f0,
     ]
-    @test gw.river.flux[1] ≈ -50.568389286368145f0
+    @test gw.river.flux[1] ≈ -51.34674583702381f0
     @test gw.drain.flux[1] ≈ 0.0
     @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
 end
@@ -93,12 +93,12 @@ model = Wflow.run_timestep(model)
 @testset "river domain (local inertial)" begin
     q = model.lateral.river.q_av
     river = model.lateral.river
-    @test sum(q) ≈ 0.026802342283209452f0
-    @test q[6] ≈ 0.005975747422547475f0
-    @test river.volume[6] ≈ 7.509678393246221f0
-    @test river.inwater[6] ≈ 0.00017989723366752878f0
-    @test q[13] ≈ 0.0004590884299495001f0
-    @test q[5] ≈ 0.006328157455390906f0
+    @test sum(q) ≈ 0.02727911500112358f0
+    @test q[6] ≈ 0.006111263175002127f0
+    @test river.volume[6] ≈ 7.6120096530771075f0
+    @test river.inwater[6] ≈ 0.00022087679662860144f0
+    @test q[13] ≈ 0.0004638698607639214f0
+    @test q[5] ≈ 0.0064668491697542786f0
 end
 Wflow.close_files(model, delete_output = false)
 
@@ -123,14 +123,14 @@ model = Wflow.run_timestep(model)
 
 @testset "river and land domain (local inertial)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 0.026809644054064f0
-    @test q[6] ≈ 0.005977572310302902f0
-    @test q[13] ≈ 0.0004591969811223254f0
-    @test q[5] ≈ 0.006330164124991123f0
+    @test sum(q) ≈ 0.027286431923384962f0
+    @test q[6] ≈ 0.00611309161099138f0
+    @test q[13] ≈ 0.0004639786629631376f0
+    @test q[5] ≈ 0.006468859889145798f0
     h = model.lateral.river.h_av
-    @test h[6] ≈ 0.08033159784710718f0
-    @test h[5] ≈ 0.0777121612660625f0
-    @test h[13] ≈ 0.08236630657766124f0
+    @test h[6] ≈ 0.08120137914886108f0
+    @test h[5] ≈ 0.07854966203902745f0
+    @test h[13] ≈ 0.08323543174453409f0
     qx = model.lateral.land.qx
     qy = model.lateral.land.qy
     @test all(qx .== 0.0f0)
@@ -152,8 +152,8 @@ model = Wflow.run_timestep(model)
 @testset "second timestep warm start" begin
     sbm = model.vertical
     @test sbm.runoff[1] == 0.0
-    @test sbm.soilevap[1] ≈ 0.2926821507238362f0
-    @test sbm.transpiration[1] ≈ 1.0122634204681036f0
+    @test sbm.soilevap[1] ≈ 0.2889306511074693f0
+    @test sbm.transpiration[1] ≈ 0.8370726722706481f0
 end
 
 @testset "overland flow warm start (kinematic wave)" begin
@@ -164,25 +164,25 @@ end
 @testset "river domain warm start (kinematic wave)" begin
     q = model.lateral.river.q_av
     river = model.lateral.river
-    @test sum(q) ≈ 0.011909898793648826f0
-    @test q[6] ≈ 0.0024332762253494573f0
-    @test river.volume[6] ≈ 2.226607059761764f0
-    @test river.inwater[6] ≈ -1.3378627340716941f-5
-    @test q[13] ≈ 7.318205070445492f-5
-    @test q[network.river.order[end]] ≈ 0.002470511395424463f0
+    @test sum(q) ≈ 0.01191742350356312f0
+    @test q[6] ≈ 0.0024353072305122064f0
+    @test river.volume[6] ≈ 2.2277585577366357f0
+    @test river.inwater[6] ≈ -1.3019072795599315f-5
+    @test q[13] ≈ 7.332742814063803f-5
+    @test q[network.river.order[end]] ≈ 0.002472526149620472f0
 end
 
 @testset "groundwater warm start" begin
     gw = model.lateral.subsurface
     @test gw.river.stage[1] ≈ 1.2031171676781156f0
     @test gw.flow.aquifer.head[17:21] ≈ [
-        1.227685769857181f0,
-        1.286975535107666f0,
+        1.2277456867225283f0,
+        1.286902494792006f0,
         1.7999999523162842f0,
-        1.5892917098077068f0,
-        1.2093841308150408f0,
+        1.5901747932190804f0,
+        1.2094238817776854f0,
     ]
-    @test gw.river.flux[1] ≈ -6.680320482009705f0
+    @test gw.river.flux[1] ≈ -6.692884222603261f0
     @test gw.drain.flux[1] ≈ 0.0
     @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
 end
diff --git a/test/run_sbm_gwf_piave.jl b/test/run_sbm_gwf_piave.jl
new file mode 100644
index 000000000..2f39b5a90
--- /dev/null
+++ b/test/run_sbm_gwf_piave.jl
@@ -0,0 +1,50 @@
+tomlpath = joinpath(@__DIR__, "sbm_gwf_piave_demand_config.toml")
+config = Wflow.Config(tomlpath)
+model = Wflow.initialize_sbm_gwf_model(config)
+model = Wflow.run_timestep(model)
+sbm = model.vertical
+
+@testset "piave water demand and allocation first timestep" begin
+    sum_total_alloc = sum(sbm.allocation.total_alloc)
+    @test sum(sbm.allocation.irri_alloc) + sum(sbm.allocation.nonirri_alloc) ≈
+          sum_total_alloc
+    @test sum(sbm.allocation.surfacewater_alloc) ≈ 1016.4268126167575f0
+    @test sum(sbm.allocation.act_groundwater_abst) ≈ 182.2057678312209f0
+    @test sbm.paddy.h[[45, 76, 296]] ≈
+          [33.55659436283413f0, 44.11663357735189f0, 35.232731550849486f0]
+    @test sbm.paddy.irrigation_trigger[[45, 76, 296]] == [1, 1, 1]
+    @test sbm.paddy.demand_gross[[45, 76, 296]] ≈ [0.0, 0.0, 0.0]
+    @test sbm.nonpaddy.irrigation_trigger[[10, 33, 1293]] == [1, 1, 1]
+    @test sbm.nonpaddy.demand_gross[[10, 33, 1293]] ≈ [3.3014913197447964f0, 0.0f0, 0.0f0]
+    @test sbm.industry.demand_gross[[1, end]] ≈ [0.2105557769536972f0, 0.0485190823674202f0]
+    @test sbm.industry.demand_net[[1, end]] ≈
+          [0.05265098437666893f0, 0.012132546864449978f0]
+    @test sbm.industry.returnflow[[1, end]] ≈ [0.15790479257702827f0, 0.03638653550297022f0]
+    @test sbm.livestock.demand_gross[[1, end]] ≈
+          [9.896758274408057f-5, 6.352497439365834f-5]
+    @test sbm.livestock.demand_net[[1, end]] ≈ [9.896758274408057f-5, 6.352497439365834f-5]
+    @test sbm.livestock.returnflow[[1, end]] ≈ [0.0f0, 0.0f0]
+    @test sbm.domestic.demand_gross[[1, end]] ≈ [0.5389957427978516f0, 0.0f0]
+    @test sbm.domestic.demand_net[[1, end]] ≈ [0.33949509263038635f0, 0.0f0]
+    @test sbm.domestic.returnflow[[1, end]] ≈ [0.1995004952035704f0, 0.0f0]
+end
+
+model = Wflow.run_timestep(model)
+sbm = model.vertical
+
+@testset "piave water demand and allocation second timestep" begin
+    sum_total_alloc = sum(sbm.allocation.total_alloc)
+    @test sum(sbm.allocation.irri_alloc) + sum(sbm.allocation.nonirri_alloc) ≈
+          sum_total_alloc
+    @test sum(sbm.allocation.surfacewater_alloc) ≈ 1057.2573971583527f0
+    @test sum(sbm.allocation.act_groundwater_abst) ≈ 189.12436534660375f0
+    @test sbm.paddy.h[[45, 76, 296]] ≈
+          [28.197082339552537f0, 25.873022895247782f0, 30.066801639786547f0]
+    @test sbm.paddy.irrigation_trigger[[45, 76, 296]] == [1, 1, 1]
+    @test sbm.paddy.demand_gross[[45, 76, 296]] ≈ [0.0f0, 0.0f0, 0.0f0]
+    @test sbm.nonpaddy.irrigation_trigger[[10, 33, 1293]] == [1, 1, 1]
+    @test sbm.nonpaddy.demand_gross[[10, 33, 1293]] ≈
+          [4.059144161330735f0, 0.0f0, 1.9399078662788196f0]
+end
+
+Wflow.close_files(model, delete_output = false)
diff --git a/test/run_sbm_piave.jl b/test/run_sbm_piave.jl
new file mode 100644
index 000000000..1d5a2fae6
--- /dev/null
+++ b/test/run_sbm_piave.jl
@@ -0,0 +1,153 @@
+
+function run_piave(model, steps)
+    q = zeros(steps)
+    ssf_vol = zeros(steps)
+    riv_vol = zeros(steps)
+    for i = 1:steps
+        model = Wflow.run_timestep(model)
+        ssf_vol[i] = mean(model.lateral.subsurface.volume)
+        riv_vol[i] = mean(model.lateral.river.volume)
+        q[i] = model.lateral.river.q_av[1]
+    end
+    return q, riv_vol, ssf_vol
+end
+
+tomlpath = joinpath(@__DIR__, "sbm_piave_demand_config.toml")
+config = Wflow.Config(tomlpath)
+model = Wflow.initialize_sbm_model(config)
+q_demand, riv_vol_demand, ssf_vol_demand = run_piave(model, 30)
+Wflow.close_files(model, delete_output = false)
+
+tomlpath = joinpath(@__DIR__, "sbm_piave_config.toml")
+config = Wflow.Config(tomlpath)
+model = Wflow.initialize_sbm_model(config)
+q_, riv_vol, ssf_vol = run_piave(model, 30)
+Wflow.close_files(model, delete_output = false)
+
+@testset "piave with and without water demand" begin
+    idx = 1:3:28
+    @test q_demand[idx] ≈ [
+        218.52770747627918f0,
+        193.02890386240665f0,
+        271.68778616960685f0,
+        161.80734173386602f0,
+        164.81279958487485f0,
+        150.4149716788231f0,
+        121.02659706677429f0,
+        108.15426854851842f0,
+        174.63022487569546f0,
+        108.20883498755789f0,
+    ]
+    @test q_[idx] ≈ [
+        219.9042687883228f0,
+        197.73238933696658f0,
+        276.99163278840734f0,
+        165.74244080863346f0,
+        172.99856395134805f0,
+        153.91963517651158f0,
+        128.52653903788593f0,
+        112.02923578000491f0,
+        178.19599851038708f0,
+        109.99414262238396f0,
+    ]
+    @test riv_vol_demand[idx] ≈ [
+        60125.790043761706f0,
+        55476.580177102805f0,
+        62195.90597660078f0,
+        48963.00368970478f0,
+        51871.46897847274f0,
+        43794.737044826295f0,
+        41710.964839545784f0,
+        39623.67826731185f0,
+        44719.74807290461f0,
+        39494.768128540185f0,
+    ]
+    @test riv_vol[idx] ≈ [
+        60620.42943610538f0,
+        55927.61042067993f0,
+        62448.21099253601f0,
+        49283.04780863544f0,
+        52293.53301649927f0,
+        44029.56815119591f0,
+        41996.885591977836f0,
+        39751.39389784692f0,
+        44619.74874948371f0,
+        39285.722489743566f0,
+    ]
+    @test ssf_vol_demand[idx] ≈ [
+        244149.81771426898f0,
+        239205.5037877743f0,
+        237749.2625556856f0,
+        233257.1175894683f0,
+        229822.63769575363f0,
+        225255.17832343685f0,
+        220944.8345415463f0,
+        216807.0855559181f0,
+        216112.96364731618f0,
+        212731.17253347262f0,
+    ]
+    @test ssf_vol[idx] ≈ [
+        244138.43308252218f0,
+        239156.33098526628f0,
+        237706.4375652105f0,
+        233210.95033208298f0,
+        229790.6284027817f0,
+        225218.1340839258f0,
+        220921.3928333043f0,
+        216784.53043204424f0,
+        216074.33950813126f0,
+        212679.70047345175f0,
+    ]
+end
+
+tomlpath = joinpath(@__DIR__, "sbm_piave_demand_config.toml")
+config = Wflow.Config(tomlpath)
+model = Wflow.initialize_sbm_model(config)
+model = Wflow.run_timestep(model)
+sbm = model.vertical
+
+@testset "piave water demand and allocation first timestep" begin
+    sum_total_alloc = sum(sbm.allocation.total_alloc)
+    @test sum(sbm.allocation.irri_alloc) + sum(sbm.allocation.nonirri_alloc) ≈
+          sum_total_alloc
+    @test sum(sbm.allocation.surfacewater_alloc) ≈ 1030.1528204311428f0
+    @test sum(sbm.allocation.act_groundwater_abst) ≈ 184.47031930837645f0
+    @test sbm.paddy.h[[45, 76, 296]] ≈
+          [34.759292485507515f0, 42.517504464353635f0, 35.83766686539591f0]
+    @test sbm.paddy.irrigation_trigger[[45, 76, 296]] == [1, 1, 1]
+    @test sbm.paddy.demand_gross[[45, 76, 296]] ≈ [0.0, 0.0, 0.0]
+    @test sbm.nonpaddy.irrigation_trigger[[10, 33, 1293]] == [1, 1, 1]
+    @test sbm.nonpaddy.demand_gross[[10, 33, 1293]] ≈
+          [3.031574420740574f0, 1.8618934872392217f0, 0.42233065562148375f0]
+    @test sbm.industry.demand_gross[[1, end]] ≈ [0.2105557769536972f0, 0.0485190823674202f0]
+    @test sbm.industry.demand_net[[1, end]] ≈
+          [0.05265098437666893f0, 0.012132546864449978f0]
+    @test sbm.industry.returnflow[[1, end]] ≈ [0.15790479257702827f0, 0.03638653550297022f0]
+    @test sbm.livestock.demand_gross[[1, end]] ≈
+          [9.896758274408057f-5, 6.352497439365834f-5]
+    @test sbm.livestock.demand_net[[1, end]] ≈ [9.896758274408057f-5, 6.352497439365834f-5]
+    @test sbm.livestock.returnflow[[1, end]] ≈ [0.0f0, 0.0f0]
+    @test sbm.domestic.demand_gross[[1, end]] ≈ [0.5389957427978516f0, 0.0f0]
+    @test sbm.domestic.demand_net[[1, end]] ≈ [0.33949509263038635f0, 0.0f0]
+    @test sbm.domestic.returnflow[[1, end]] ≈ [0.1995004952035704f0, 0.0f0]
+end
+
+model = Wflow.run_timestep(model)
+sbm = model.vertical
+
+@testset "piave water demand and allocation second timestep" begin
+    sum_total_alloc = sum(sbm.allocation.total_alloc)
+    @test sum(sbm.allocation.irri_alloc) + sum(sbm.allocation.nonirri_alloc) ≈
+          sum_total_alloc
+    @test sum(sbm.allocation.surfacewater_alloc) ≈ 940.1941924010235f0
+    @test sum(sbm.allocation.act_groundwater_abst) ≈ 163.59123515939052f0
+    @test sbm.paddy.h[[45, 76, 296]] ≈
+          [29.50674105890283f0, 38.11966817469463f0, 30.679920457042897f0]
+    @test sbm.paddy.irrigation_trigger[[45, 76, 296]] == [1, 1, 1]
+    @test sbm.paddy.demand_gross[[45, 76, 296]] ≈ [0.0f0, 0.0f0, 0.0f0]
+    @test sbm.nonpaddy.irrigation_trigger[[10, 33, 1293]] == [1, 1, 1]
+    @test sbm.nonpaddy.demand_gross[[10, 33, 1293]] ≈
+          [3.8551909476218054f0, 0.0f0, 1.3531153828385536f0]
+end
+
+Wflow.close_files(model, delete_output = false)
diff --git a/test/runtests.jl b/test/runtests.jl
index 18771d1de..3dd096afe 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -70,7 +70,10 @@ lake_sh_2_path = testdata(v"0.2.1", "lake_sh_2.csv", "lake_sh_2.csv")
 lake_hq_2_path = testdata(v"0.2.1", "lake_hq_2.csv", "lake_hq_2.csv")
 forcing_calendar_noleap_path =
     testdata(v"0.2.8", "forcing-calendar-noleap.nc", "forcing-calendar-noleap.nc")
-
+forcing_piave_path = testdata(v"0.2.9", "inmaps-era5-2010-piave.nc", "forcing-piave.nc")
+staticmaps_piave_path = testdata(v"0.2.9", "staticmaps-piave.nc", "staticmaps-piave.nc")
+instates_piave_path = testdata(v"0.2.9", "instates-piave.nc", "instates-piave.nc")
+instates_piave_gwf_path = testdata(v"0.2.9", "instates-piave-gwf.nc", "instates-piave-gwf.nc")
 
 include("testing_utils.jl")
 
@@ -85,6 +88,8 @@ with_logger(NullLogger()) do
         include("vertical_process.jl")
         include("reservoir_lake.jl")
         include("run_sbm.jl")
+        include("run_sbm_piave.jl")
+        include("run_sbm_gwf_piave.jl")
         include("run_hbv.jl")
         include("run_sbm_gwf.jl")
         include("run.jl")
diff --git a/test/sbm_gwf_piave_demand_config.toml b/test/sbm_gwf_piave_demand_config.toml
new file mode 100644
index 000000000..fa7e4cb44
--- /dev/null
+++ b/test/sbm_gwf_piave_demand_config.toml
@@ -0,0 +1,189 @@
+calendar = "proleptic_gregorian"
+starttime = 2010-07-01T00:00:00
+endtime = 2010-10-01T00:00:00
+time_units = "days since 1900-01-01 00:00:00"
+timestepsecs = 86400
+dir_input = "data/input"
+dir_output = "data/output"
+loglevel = "info"
+
+[state]
+path_input = "instates-piave-gwf.nc"
+path_output = "outstates-piave-gwf.nc"
+
+[input]
+path_forcing = "forcing-piave.nc"
+path_static = "staticmaps-piave.nc"
+ldd = "wflow_ldd"
+river_location = "wflow_river"
+altitude = "wflow_dem"
+subcatchment = "wflow_subcatch"
+forcing = [ "vertical.precipitation", "vertical.temperature", "vertical.potential_evaporation",]
+cyclic = [ "vertical.leaf_area_index", "vertical.domestic.demand_gross", "vertical.domestic.demand_net", "vertical.industry.demand_gross", "vertical.industry.demand_net", "vertical.livestock.demand_gross", "vertical.livestock.demand_net", "vertical.paddy.irrigation_trigger", "vertical.nonpaddy.irrigation_trigger",]
+gauges = "wflow_gauges"
+gauges_grdc = "wflow_gauges_grdc"
+
+[model]
+type = "sbm_gwf"
+constanthead = true
+masswasting = true
+snow = true
+reinit = false
+reservoirs = false
+lakes = false
+glacier = true
+kin_wave_iteration = true
+kw_river_tstep = 900
+kw_land_tstep = 3600
+thicknesslayers = [ 50, 100, 50, 200, 800,]
+river_routing = "kinematic-wave"
+
+[state.vertical]
+satwaterdepth = "satwaterdepth"
+snow = "snow"
+tsoil = "tsoil"
+ustorelayerdepth = "ustorelayerdepth"
+snowwater = "snowwater"
+canopystorage = "canopystorage"
+glacierstore = "glacierstore"
+
+[state.vertical.paddy]
+h = "h_paddy"
+
+[state.lateral.subsurface.flow.aquifer]
+head = "head"
+
+[input.vertical]
+alpha_h1 = "alpha_h1"
+altitude = "wflow_dem"
+c = "c"
+cf_soil = "cf_soil"
+cfmax = "Cfmax"
+e_r = "EoverR"
+f = "f"
+infiltcappath = "InfiltCapPath"
+infiltcapsoil = "InfiltCapSoil"
+kext = "Kext"
+kv_0 = "KsatVer"
+leaf_area_index = "LAI"
+m = "M_"
+maxleakage = "MaxLeakage"
+pathfrac = "PathFrac"
+potential_evaporation = "pet"
+precipitation = "precip"
+rootdistpar = "rootdistpar"
+rootingdepth = "RootingDepth"
+soilminthickness = "SoilMinThickness"
+soilthickness = "SoilThickness_gw"
+specific_leaf = "Sl"
+storage_wood = "Swood"
+temperature = "temp"
+tt = "TT"
+tti = "TTI"
+ttm = "TTM"
+water_holding_capacity = "WHC"
+waterfrac = "WaterFrac"
+theta_s = "thetaS"
+theta_r = "thetaR"
+glacierstore = "wflow_glacierstore"
+glacierfrac = "wflow_glacierfrac"
+g_cfmax = "G_Cfmax"
+g_tt = "G_TT"
+g_sifrac = "G_SIfrac"
+kc = "crop_factor"
+kvfrac = "kvfrac"
+h1 = "h1"
+h2 = "h2"
+h3_high = "h3_high"
+h3_low = "h3_low"
+h4 = "h4"
+
+[model.water_demand]
+domestic = true
+industry = true
+livestock = true
+paddy = true
+nonpaddy = true
+
+[state.lateral.river]
+q = "q_river"
+h = "h_river"
+h_av = "h_av_river"
+
+[state.lateral.subsurface]
+ssf = "ssf"
+
+[state.lateral.land]
+q = "q_land"
+h = "h_land"
+h_av = "h_av_land"
+
+[input.vertical.allocation]
+areas = "allocation_areas"
+frac_sw_used = "SurfaceWaterFrac"
+
+[input.vertical.domestic]
+demand_gross = "dom_gross"
+demand_net = "dom_net"
+
+[input.vertical.industry]
+demand_gross = "ind_gross"
+demand_net = "ind_net"
+
+[input.vertical.livestock]
+demand_gross = "lsk_gross"
+demand_net = "lsk_net"
+
+[input.vertical.paddy]
+irrigation_areas = "paddy_irrigation_areas"
+irrigation_trigger = "irrigation_trigger"
+
+[input.vertical.nonpaddy]
+irrigation_areas = "nonpaddy_irrigation_areas"
+irrigation_trigger = "irrigation_trigger"
+
+[input.lateral.river]
+length = "wflow_riverlength"
+n = "N_River"
+slope = "RiverSlope"
+width = "wflow_riverwidth"
+bankfull_depth = "RiverDepth"
+
+[input.lateral.subsurface]
+constant_head = "constant_head"
+conductivity_profile = "exponential"
+conductivity = "kh_surface"
+exfiltration_conductance = "riverbed_cond"
+infiltration_conductance = "riverbed_cond"
+river_bottom = "zb_river"
+specific_yield = "specific_yield"
+gwf_f = "gwf_f"
+
+[input.lateral.land]
+n = "N"
+slope = "Slope"
+
+[output]
+path = "output-piave-gwf.nc"
+
+[output.lateral.river]
+q_av = "q_river"
+
+[output.vertical]
+zi = "zi"
+
+[output.lateral.subsurface.flow.aquifer]
+head = "head"
+
+[csv]
+path = "output.csv"
+
+[[csv.column]]
+header = "Q"
+map = "gauges"
+parameter = "lateral.river.q_av"
+
+[[csv.column]]
+header = "Q"
+map = "gauges_grdc"
+parameter = "lateral.river.q_av"
\ No newline at end of file
diff --git a/test/sbm_piave_config.toml b/test/sbm_piave_config.toml
new file mode 100644
index 000000000..22c5e30c8
--- /dev/null
+++ b/test/sbm_piave_config.toml
@@ -0,0 +1,136 @@
+calendar = "proleptic_gregorian"
+starttime = 2010-07-01T00:00:00
+endtime = 2010-10-01T00:00:00
+time_units = "days since 1900-01-01 00:00:00"
+timestepsecs = 86400
+dir_input = "data/input"
+dir_output = "data/output"
+loglevel = "info"
+
+[state]
+path_input = "instates-piave.nc"
+path_output = "outstates-piave.nc"
+
+[input]
+path_forcing = "forcing-piave.nc"
+path_static = "staticmaps-piave.nc"
+ldd = "wflow_ldd"
+river_location = "wflow_river"
+subcatchment = "wflow_subcatch"
+forcing = [ "vertical.precipitation", "vertical.temperature", "vertical.potential_evaporation",]
+cyclic = [ "vertical.leaf_area_index",]
+gauges = "wflow_gauges"
+gauges_grdc = "wflow_gauges_grdc"
+
+[model]
+type = "sbm"
+masswasting = true
+snow = true
+reinit = false
+reservoirs = false
+lakes = false
+glacier = true
+kin_wave_iteration = true
+kw_river_tstep = 900
+kw_land_tstep = 3600
+thicknesslayers = [ 50, 100, 50, 200, 800,]
+river_routing = "kinematic-wave"
+
+[state.vertical]
+satwaterdepth = "satwaterdepth"
+snow = "snow"
+tsoil = "tsoil"
+ustorelayerdepth = "ustorelayerdepth"
+snowwater = "snowwater"
+canopystorage = "canopystorage"
+glacierstore = "glacierstore"
+
+[input.vertical]
+alpha_h1 = "alpha_h1"
+altitude = "wflow_dem"
+c = "c"
+cf_soil = "cf_soil"
+cfmax = "Cfmax"
+e_r = "EoverR"
+f = "f"
+infiltcappath = "InfiltCapPath"
+infiltcapsoil = "InfiltCapSoil"
+kext = "Kext"
+kv_0 = "KsatVer"
+leaf_area_index = "LAI"
+m = "M_"
+maxleakage = "MaxLeakage"
+pathfrac = "PathFrac"
+potential_evaporation = "pet"
+precipitation = "precip"
+rootdistpar = "rootdistpar"
+rootingdepth = "RootingDepth"
+soilminthickness = "SoilMinThickness"
+soilthickness = "SoilThickness"
+specific_leaf = "Sl"
+storage_wood = "Swood"
+temperature = "temp"
+tt = "TT"
+tti = "TTI"
+ttm = "TTM"
+water_holding_capacity = "WHC"
+waterfrac = "WaterFrac"
+theta_s = "thetaS"
+theta_r = "thetaR"
+glacierstore = "wflow_glacierstore"
+glacierfrac = "wflow_glacierfrac"
+g_cfmax = "G_Cfmax"
+g_tt = "G_TT"
+g_sifrac = "G_SIfrac"
+kc = "crop_factor"
+h1 = "h1"
+h2 = "h2"
+h3_high = "h3_high"
+h3_low = "h3_low"
+h4 = "h4"
+
+[state.lateral.river]
+q = "q_river"
+h = "h_river"
+h_av = "h_av_river"
+
+[state.lateral.subsurface]
+ssf = "ssf"
+
+[state.lateral.land]
+q = "q_land"
+h = "h_land"
+h_av = "h_av_land"
+
+[input.lateral.river]
+length = "wflow_riverlength"
+n = "N_River"
+slope = "RiverSlope"
+width = "wflow_riverwidth"
+bankfull_depth = "RiverDepth"
+
+[input.lateral.subsurface]
+ksathorfrac = "KsatHorFrac"
+
+[input.lateral.land]
+n = "N"
+slope = "Slope"
+
+[output]
+path = "output-piave.nc"
+
+[output.lateral.river]
+q_av = "q_river"
+
+[csv]
+path = "output-piave.csv"
+
+[[csv.column]]
+header = "Q"
+map = "gauges"
+parameter = "lateral.river.q_av"
+
+[[csv.column]]
+header = "Q"
+map = "gauges_grdc"
+parameter = "lateral.river.q_av"
\ No newline at end of file
diff --git a/test/sbm_piave_demand_config.toml b/test/sbm_piave_demand_config.toml
new file mode 100644
index 000000000..141ef0105
--- /dev/null
+++ b/test/sbm_piave_demand_config.toml
@@ -0,0 +1,186 @@
+calendar = "proleptic_gregorian"
+starttime = 2010-07-01T00:00:00
+endtime = 2010-10-01T00:00:00
+time_units = "days since 1900-01-01 00:00:00"
+timestepsecs = 86400
+dir_input = "data/input"
+dir_output = "data/output"
+loglevel = "info"
+
+[state]
+path_input = "instates-piave.nc"
+path_output = "outstates-piave.nc"
+
+[input]
+path_forcing = "forcing-piave.nc"
+path_static = "staticmaps-piave.nc"
+ldd = "wflow_ldd"
+river_location = "wflow_river"
+subcatchment = "wflow_subcatch"
+forcing = [ "vertical.precipitation", "vertical.temperature", "vertical.potential_evaporation",]
+cyclic = [ "vertical.leaf_area_index", "vertical.domestic.demand_gross", "vertical.domestic.demand_net", "vertical.industry.demand_gross", "vertical.industry.demand_net", "vertical.livestock.demand_gross", "vertical.livestock.demand_net", "vertical.paddy.irrigation_trigger", "vertical.nonpaddy.irrigation_trigger",]
+gauges = "wflow_gauges"
+gauges_grdc = "wflow_gauges_grdc"
+
+[model]
+type = "sbm"
+masswasting = true
+snow = true
+reinit = false
+reservoirs = false
+lakes = false
+glacier = true
+kin_wave_iteration = true
+kw_river_tstep = 900
+kw_land_tstep = 3600
+thicknesslayers = [ 50, 100, 50, 200, 800,]
+river_routing = "kinematic-wave"
+
+[state.vertical]
+satwaterdepth = "satwaterdepth"
+snow = "snow"
+tsoil = "tsoil"
+ustorelayerdepth = "ustorelayerdepth"
+snowwater = "snowwater"
+canopystorage = "canopystorage"
+glacierstore = "glacierstore"
+
+[state.vertical.paddy]
+h = "h_paddy"
+
+[input.vertical]
+alpha_h1 = "alpha_h1"
+altitude = "wflow_dem"
+c = "c"
+cf_soil = "cf_soil"
+cfmax = "Cfmax"
+e_r = "EoverR"
+f = "f"
+infiltcappath = "InfiltCapPath"
+infiltcapsoil = "InfiltCapSoil"
+kext = "Kext"
+kv_0 = "KsatVer"
+leaf_area_index = "LAI"
+m = "M_"
+maxleakage = "MaxLeakage"
+pathfrac = "PathFrac"
+potential_evaporation = "pet"
+precipitation = "precip"
+rootdistpar = "rootdistpar"
+rootingdepth = "RootingDepth"
+soilminthickness = "SoilMinThickness"
+soilthickness = "SoilThickness"
+specific_leaf = "Sl"
+storage_wood = "Swood"
+temperature = "temp"
+tt = "TT"
+tti = "TTI"
+ttm = "TTM"
+water_holding_capacity = "WHC"
+waterfrac = "WaterFrac"
+theta_s = "thetaS"
+theta_r = "thetaR"
+glacierstore = "wflow_glacierstore"
+glacierfrac = "wflow_glacierfrac"
+g_cfmax = "G_Cfmax"
+g_tt = "G_TT"
+g_sifrac = "G_SIfrac"
+kc = "crop_factor"
+kvfrac = "kvfrac"
+h1 = "h1"
+h2 = "h2"
+h3_high = "h3_high"
+h3_low = "h3_low"
+h4 = "h4"
+
+[model.water_demand]
+domestic = true
+industry = true
+livestock = true
+paddy = true
+nonpaddy = true
+
+[state.lateral.river]
+q = "q_river"
+h = "h_river"
+h_av = "h_av_river"
+
+[state.lateral.subsurface]
+ssf = "ssf"
+
+[state.lateral.land]
+q = "q_land"
+h = "h_land"
+h_av = "h_av_land"
+
+[input.vertical.allocation]
+areas = "allocation_areas"
+frac_sw_used = "SurfaceWaterFrac"
+
+[input.vertical.domestic]
+demand_gross = "dom_gross"
+demand_net = "dom_net"
+
+[input.vertical.industry]
+demand_gross = "ind_gross"
+demand_net = "ind_net"
+
+[input.vertical.livestock]
+demand_gross = "lsk_gross"
+demand_net = "lsk_net"
+
+[input.vertical.paddy]
+irrigation_areas = "paddy_irrigation_areas"
+irrigation_trigger = "irrigation_trigger"
+
+[input.vertical.nonpaddy]
+irrigation_areas = "nonpaddy_irrigation_areas"
+irrigation_trigger = "irrigation_trigger"
+
+[input.lateral.river]
+length = "wflow_riverlength"
+n = "N_River"
+slope = "RiverSlope"
+width = "wflow_riverwidth"
+bankfull_depth = "RiverDepth"
+
+[input.lateral.subsurface]
+ksathorfrac = "KsatHorFrac"
+
+[input.lateral.land]
+n = "N"
+slope = "Slope"
+
+[output]
+path = "output-piave-demand.nc"
+
+[output.lateral.river]
+q_av = "q_river"
+
+[output.vertical]
+zi = "zi"
+
+[csv]
+path = "output-piave-demand.csv"
+
+[[csv.column]]
+header = "Q"
+map = "gauges"
+parameter = "lateral.river.q_av"
+
+[[csv.column]]
+header = "Q"
+map = "gauges_grdc"
+parameter = "lateral.river.q_av"
+
+[[csv.column]]
+coordinate.x = 12.7243
+coordinate.y = 45.5851
+header = "paddy_h_bycoord"
+parameter = "vertical.paddy.h"
+
+[[csv.column]]
+coordinate.x = 12.7243
+coordinate.y = 45.5851
+header = "irri_bycoord"
+parameter = "vertical.allocation.irri_alloc"
\ No newline at end of file
diff --git a/test/utils.jl b/test/utils.jl
index 3edf84c19..c1922fb6d 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -21,3 +21,11 @@ end
     @test Wflow.julian_day(DateTime(2001, 3, 1)) == 60
     @test Wflow.julian_day(DateTime(2001, 12, 31)) == 365
 end
+
+@testset "Bounded divide" begin
+    @test Wflow.bounded_divide(1.0, 0.0) == 0.0
+    @test Wflow.bounded_divide(1.0, 0.0, default = 0.5) == 0.5
+    @test Wflow.bounded_divide(1.0, 0.5) == 1.0
+    @test Wflow.bounded_divide(1.0, 0.5, max = 0.75) == 0.75
+    @test Wflow.bounded_divide(1.0, 2.0) == 0.5
+end
diff --git a/test/vertical_process.jl b/test/vertical_process.jl
index 7118a8b24..9ce59fb45 100644
--- a/test/vertical_process.jl
+++ b/test/vertical_process.jl
@@ -1,3 +1,4 @@
+using Dates
 @testset "vertical processes" begin
     @test all(
         isapprox.(
@@ -11,28 +12,22 @@
             (4.343, 3.87, 0.387, 0.0, 3.8, 2.043),
         ),
     )
-    @test all(
-        isapprox.(
-            Wflow.acttransp_unsat_sbm(
-                300.0,
-                55.0,
-                400.0,
-                3.6,
-                1.85,
-                10.5,
-                300.0,
-                0.6,
-                0.1,
-                10.0,
-                false,
-            ),
-            (55.0, 1.85, 3.6),
-        ),
-    )
+    @test Wflow.head_brooks_corey(0.25, 0.6, 0.15, 10.5, -10.0) ≈ -90.6299820833844
+    @test Wflow.feddes_h3(-300.0, -600.0, 3.5, Second(86400)) ≈ -412.5
+    @test Wflow.feddes_h3(-300.0, -600.0, 0.5, Second(86400)) == -600.0
+    @test Wflow.feddes_h3(-300.0, -600.0, 6.0, Second(86400)) == -300.0
+    @test Wflow.rwu_reduction_feddes(0.0, -10.0, -100.0, -300.0, -15000.0, 0.0) == 0.0
+    @test Wflow.rwu_reduction_feddes(0.0, -10.0, -100.0, -300.0, -15000.0, 1.0) == 1.0
+    @test Wflow.rwu_reduction_feddes(-90.0, -10.0, -100.0, -412.5, -15000.0, 0.0) ≈
+          0.8888888888888888
+    @test Wflow.rwu_reduction_feddes(-350.0, -10.0, -100.0, -412.5, -15000.0, 0.0) == 1.0
+    @test Wflow.rwu_reduction_feddes(-12000.0, -10.0, -100.0, -412.5, -15000.0, 0.0) ≈
+          0.20565552699228792
+    @test Wflow.rwu_reduction_feddes(-16000.0, -10.0, -100.0, -412.5, -15000.0, 0.0) == 0.0
     @test all(
         isapprox.(
             Wflow.infiltration(27.5, 0.2, 0.038, 8.9, 50.0, 5.0, 23.5, false, false),
-            (23.5, 19.14814814814815, 4.351851851851852, 22.0, 5.5, 0.5),
+            (23.5, 19.14814814814815, 4.351851851851852, 22.0, 5.5, 0.5, 1.0),
         ),
     )
     @test all(