diff --git a/src/Model/Geometry/Profile.py b/src/Model/Geometry/Profile.py
index ad1b82252ab790872796e8f402ac428aea0ef438..ee662308f33def332297ed0d450d552d74b30669 100644
--- a/src/Model/Geometry/Profile.py
+++ b/src/Model/Geometry/Profile.py
@@ -351,3 +351,7 @@ class Profile(object):
     # Abstract method for width approximation
     def width_approximation(self):
         raise NotImplementedMethodeError(self, self.width_approximation)
+
+    # Abstract method for width approximation
+    def get_water_limits(self, z):
+        raise NotImplementedMethodeError(self, self.get_water_limits)
diff --git a/src/Model/Geometry/ProfileXYZ.py b/src/Model/Geometry/ProfileXYZ.py
index dfe187024d3c8cee425c54adfa07e09985a47152..2763924b72a1ce561347defcb161fddacda2aa77 100644
--- a/src/Model/Geometry/ProfileXYZ.py
+++ b/src/Model/Geometry/ProfileXYZ.py
@@ -442,3 +442,44 @@ class ProfileXYZ(Profile, SQLSubModel):
         rd = self.get_point_by_name("rd")
 
         return abs(rg.dist(rd))
+
+    def get_water_limits(self, z):
+#==============================================================================
+#   détermination des points limites RG et RD pour un niveau d'eau donné
+#
+#   irg et ird sont les premiers indices en partant des rives gauche et
+#   droite qui correspondent à des points sous la surface de l'eau
+#   ptX et ptY sont les points interpolés où le plan d'eau intersecte le profil
+#   known_level est le niveau d'eau pour lequel on a obtenu irg, ird, ptX et ptY
+#==============================================================================
+        # initialisation
+        irg = -1  ;  ird = -1
+        for i in range(self.number_points):
+            if self.point(i).z <= z:
+                irg = i
+        for i in reversed(range(self.number_points)):
+            if self.point(i).z <= z:
+                ird = i
+        # interpolation des points ptX et ptY
+        if (irg > 0):
+            x=np.interp(z,
+                        [self.point(irg).z,self.point(irg+1).z],
+                        [self.point(irg).x,self.point(irg+1).x])
+            y=np.interp(z,
+                        [self.point(irg).z,self.point(irg+1).z],
+                        [self.point(irg).y,self.point(irg+1).y])
+            ptX=PointXYZ(x,y,z)
+        else:
+            ptX = self.point(0)
+        if (ird < self.number_points-1):
+            x=np.interp(z,
+                        [self.point(ird).z,self.point(ird+1).z],
+                        [self.point(ird).x,self.point(ird+1).x])
+            y=np.interp(z,
+                        [self.point(ird).z,self.point(ird+1).z],
+                        [self.point(ird).y,self.point(ird+1).y])
+            ptY=PointXYZ(x,y,z)
+        else:
+            ptY = self.point(self.number_points-1)
+
+        return ptX,ptY
diff --git a/src/Solver/Mage.py b/src/Solver/Mage.py
index fa4ce18c1369d0f539c57c37b97cb0e266115060..22349870d0e3da5b1ae7f644ff4e984d551deae1 100644
--- a/src/Solver/Mage.py
+++ b/src/Solver/Mage.py
@@ -813,6 +813,12 @@ class Mage8(Mage):
 
                         # Set data for profile RI
                         reach.set(ri, timestamp, key, d)
+                        if key == "Z":
+                            profile = study.river.current_reach().reach.profile(ri)
+                            ptX,ptY = profile.get_water_limits(d)
+                            reach.set(ri, timestamp, "ptX", ptX)
+                            reach.set(ri, timestamp, "ptY", ptY)
+
 
                 endline()
                 end = newline().size <= 0