Index: wflow-py/wflow/wflow_sbm.py
===================================================================
diff -u -rf33154035afcaa101c5a30594b9c86a4fe9f61ef -re3d3271bfc58a9021d816e3332479f5a040afdc8
--- wflow-py/wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision f33154035afcaa101c5a30594b9c86a4fe9f61ef)
+++ wflow-py/wflow/wflow_sbm.py	(.../wflow_sbm.py)	(revision e3d3271bfc58a9021d816e3332479f5a040afdc8)
@@ -2892,22 +2892,43 @@
             Runoff = self.SurfaceRunoff
 
         # Determine Soil moisture profile
-        # 1: average volumetric soil in total unsat store
-        # self.SMVol = (cover(self.UStoreDepth/self.zi,0.0) + self.thetaR) * (self. thetaS - self.thetaR)
-        # self.SMRootVol = (cover(self.UStoreDepth/min(self.ActRootingDepth,self.zi),0.0) + self.thetaR) * (self. thetaS - self.thetaR)
-        RootStore_sat = max(0.0, self.ActRootingDepth - self.zi) * self.thetaS
-        RootStore_unsat = self.ZeroMap
+        # self.vwc, self.vwcRoot: volumetric water content [m3/m3] per soil layer and root zone (including thetaR and saturated store)
+        # self.vwc_perc, self.vwc_percRoot: volumetric water content [%] per soil layer and root zone (including thetaR and saturated store)
+        # self.RootStore_sat: root water storage [mm] in saturated store (excluding thetaR)
+        # self.RootStore_unsat: root water storage [mm] in unsaturated store (excluding thetaR)
+        # self.RootStore: total root water storage [mm] (excluding thetaR)
 
+        self.RootStore_sat = max(0.0, self.ActRootingDepth - self.zi) * (
+            self.thetaS - self.thetaR
+        )
+
         self.RootStore_unsat = self.ZeroMap
         self.SumThickness = self.ZeroMap
         self.vwc = []
         self.vwc_perc = []
 
         for n in arange(len(self.UStoreLayerThickness)):
-            RootStore_unsat = RootStore_unsat + (
-                max(0.0, (self.ActRootingDepth - self.SumThickness))
-                / (self.UStoreLayerThickness[n])
-            ) * (cover(self.UStoreLayerDepth[n] + self.thetaR, 0.0))
+
+            fracRoot = ifthenelse(
+                self.ZiLayer > n,
+                min(
+                    1.0,
+                    max(
+                        0.0,
+                        (min(self.ActRootingDepth, self.zi) - self.SumThickness)
+                        / self.UStoreLayerThickness[n],
+                    ),
+                ),
+                min(
+                    1.0,
+                    max(
+                        0.0,
+                        (self.ActRootingDepth - self.SumThickness)
+                        / (self.zi + 1 - self.SumThickness),
+                    ),
+                ),
+            )
+
             self.SumThickness = self.UStoreLayerThickness[n] + self.SumThickness
 
             self.vwc.append(