diff --git a/.gitignore b/.gitignore
index 91cc5cb2..380068f7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -77,4 +77,6 @@ doc/_build
 
 # baugetfa
 tmp*
+test/data/test_rsml_io.rsml
+out.csv
 
diff --git a/doc/example/example_rsml_io.py b/doc/example/example_rsml_io.py
index 0677e800..b86879b4 100644
--- a/doc/example/example_rsml_io.py
+++ b/doc/example/example_rsml_io.py
@@ -13,11 +13,11 @@
   into continuous mtg representation
 - set the resolution according to the unit precised in the rsml, in Hydroroot lengths are in meter
 - transform the continuous mtg to discrete mtg usable in Hydroroot according to resolution
-- calculation of g properties (radius, mylength) and  flux
+- calculation of g properties (radius, dist_to_base) and  flux
 
 - export discrete mtg to rsml
 - re-import rsml to continuous mtg, transform the latter to discrete mtg
-- re-do calculation of g properties (radius, mylength) and  flux
+- re-do calculation of g properties (radius, dist_to_base) and  flux
 
 - compare both calculations should be zero
 
@@ -53,7 +53,7 @@ def set_mtg_properties(g):
     Set MTG properties and perform some gemetrical calculation
 
     The vertex radius properties is set.
-    The following properties are computed: length, position, mylength, surface, volume, total length,
+    The following properties are computed: length, position, dist_to_base, surface, volume, total length,
         primary root length
 
     :param:
@@ -75,11 +75,11 @@ def set_mtg_properties(g):
     # Calculation of the distance from base of each vertex, used for cut and flow
     # Remark: this calculation is done in flux.segments_at_length; analysis.nb_roots but there is a concern with the
     # parameter dl which should be equal to vertex length but which is not pass
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + parameter.archi['segment_length'] if pid else parameter.archi['segment_length']
-    g.properties()['mylength'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + parameter.archi['segment_length'] if pid else parameter.archi['segment_length']
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * parameter.archi['segment_length']
@@ -139,7 +139,7 @@ def hydro_calculation(g, axfold = 1., radfold = 1., axial_data = None, k_radial
     # continuous mtg to discrete mtg
     g = import_rsml_to_discrete_mtg(g_c, segment_length = parameter.archi['segment_length'], resolution = resolution)
 
-    # calculation of g properties: radius, mylength, etc.
+    # calculation of g properties: radius, dist_to_base, etc.
     g, primary_length, _length, surface = set_mtg_properties(g)
 
     # flux calculation
diff --git a/doc/example/metafspm/example_metafspm_archi_from_file.py b/doc/example/metafspm/example_metafspm_archi_from_file.py
new file mode 100644
index 00000000..9216abe3
--- /dev/null
+++ b/doc/example/metafspm/example_metafspm_archi_from_file.py
@@ -0,0 +1,166 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# # Measured architecture
+
+# Read architecture from text file and run a simulation of the sap flux from an Arabidopsis de-topped root plunged in a hydroponic solution at a hydrostatic pressure of 0.4 Mpa when its  base is at the atmospheric pressure.
+
+# ## The text format architecture file
+# 
+# |distance_from_base_(mm)  |lateral_root_length_(mm)  |order|
+# |:-:  | :-:  | :-: |
+# |0.89                     |90.81             	       |1|
+# |3.02                     |63.98             	       |1|
+# |102.94                   |0.0             	         |1|
+# |2.14                     |23.72             	       |1-1|
+# |90.81                    |0.0             	         |1-1|
+# |2.48                     |5.15             	       |1-2|
+# |63.98                    |0.0             	         |1-2|
+# 
+# This is a tab separated text file with 3 columns:
+# 
+#     1. the distance from base of the branching laterals in mm
+# 
+#     2. the lateral root length in mm
+# 
+#     3. a string of one or more number indicating the parent root
+# 
+# In the example above, the root has two lateral of 1st order and on each of them one lateral of 2d order. The order of 1 indicates that the laterals are on the primary root. The last line with order 1 with 0.0 in the second column indicates the primary root tip. The line with order 1-1 indicates that this a second order lateral on the first lateral positioned at 2.14 mm from the branching on the primary root. And so on.
+# 
+# A 4th column with the averaged diameter of the root may be given, that may be used to build the MTG representing the architecture.
+
+# ## Running the calculation
+# If the package HydroRoot is not installed, the following examples can be run by cloning the sources from git and then sourcing the src directory in Ipython console for instance like this:
+
+# In[1]:
+
+
+get_ipython().run_line_magic('matplotlib', 'inline')
+
+
+# In[2]:
+
+
+import matplotlib
+from openalea.widgets.plantgl import PlantGL # notebook viewer 3D
+from openalea.plantgl.algo.view import view # 2D view
+from openalea.hydroroot.display import mtg_scene
+from openalea.hydroroot.read_file import read_archi_data
+from openalea.hydroroot.main import hydroroot_flow, root_builder
+from openalea.hydroroot.generator import measured_root
+from openalea.hydroroot import model
+from openalea.mtg import MTG
+
+
+# Read the architecture file as DataFrame
+
+# In[3]:
+
+
+df = read_archi_data('../data/plant-01.txt')
+
+
+# Building the MTG from the file, and  return the primary root length, the total length and the surface. The seed refer to the seed of the root generator when the MTG is not built from a file but is generated.
+
+# In[4]:
+
+segment_length = 1.0e-4
+time_step = 1.0e-4
+
+g = measured_root.mtg_from_aqua_data(df, segment_length)
+hydromodel = model.HydroRootModel(g, time_step)
+hydromodel.compute_metrics()
+
+# Some conductance data versus distance to tip
+
+# In[5]:
+
+
+hydromodel.compute_radial_conductance(data = ([0, 0.2],[30.0,30.0]))
+hydromodel.compute_axial_conductance(data = ([0, 0.2],[3.0e-7,4.0e-4]))
+
+
+# Flux and equivalent conductance calculation, for a root in an external hydroponic medium at 0.4 MPa, its base at 0.1 MPa, and with the conductances set above.
+
+# In[6]:
+
+# g = hydromodel.g
+# g, keq, jv = hydroroot_flow(g, psi_e = 0.4, psi_base = 0.1, axial_conductivity_data = K_axial_data, radial_conductivity_data = k_radial_data)
+hydromodel.psi_e = 0.4
+hydromodel.psi_base = 0.1
+hydromodel.hydrostatic_solver_flux()
+
+# In[7]:
+
+
+print('equivalent root conductance (microL/s/MPa): ',hydromodel.keq, 'sap flux (microL/s): ', hydromodel.Jv)
+
+
+# ## Plots
+
+# Display the local water uptake heatmap in 3D
+
+# In[8]:
+
+
+s = mtg_scene(g, prop_cmap = 'j') # create a scene from the mtg with the property j is the radial flux in ul/s
+PlantGL(s) # display the root into the plantgl oawidget
+
+
+# You may change the property to display to the hydrostatic pressure inside the xylem vessels for instance
+# 
+# to reduce notebook size we use here a 2D view but you can use the openalea.widgets `PlantGL(s)` to display an interactive 3D view
+
+# In[9]:
+
+
+s = mtg_scene(g, prop_cmap='psi_in')
+view(s) #PlantGL(s)
+
+
+# You may change the radial conductivity and see the impact on the water uptake
+
+# In[10]:
+
+
+hydromodel.compute_radial_conductance(data = ([0, 0.2],[300.0,300.0]))
+hydromodel.hydrostatic_solver_flux()
+print('sap flux (microL/s): ', hydromodel.Jv)
+s = mtg_scene(g, prop_cmap='j')
+view(s) # to reduce notebook size we use here a 2D view but use PlantGL(s) to 3D
+
+
+# Or the axial conductance
+
+# In[11]:
+
+
+hydromodel.compute_radial_conductance(data = ([0, 0.2],[30.0,30.0]))
+hydromodel.compute_axial_conductance(data = ([0, 0.2],[3.0e-7,4.0e-4]))
+hydromodel.hydrostatic_solver_flux()
+print('sap flux (microL/s): ', hydromodel.Jv)
+s = mtg_scene(g, prop_cmap='j')
+view(s) # to reduce notebook size we use here a 2D view but use PlantGL(s) to 3D
+
+
+# ## PlantGl viewer
+# 
+# This is only working on **local machine**.
+# 
+# HydroRoot has also a function using the standalone Plantgl viewer allowing to specify the mtg property.
+# 
+# The following will open a standalone window with a 3D representation of the root. Uncomment the lines to run the cell locally.
+
+# In[12]:
+
+
+# from openalea.hydroroot.display import plot
+# %gui qt
+# plot(g, prop_cmap='psi_in')
+from openalea.hydroroot.init_parameter import Parameters
+parameter = Parameters()
+parameter.read_file('parameters_plant_01.yml')
+scenario = parameter.metafspm_scenario()
+hydromodel = model.HydroRootModel(g, time_step, **scenario)
+hydromodel.hydrostatic_solver_flux()
+print('sap flux (microL/s): ', hydromodel.Jv)
diff --git a/doc/example/metafspm/parameters_plant_01.yml b/doc/example/metafspm/parameters_plant_01.yml
new file mode 100644
index 00000000..815f7326
--- /dev/null
+++ b/doc/example/metafspm/parameters_plant_01.yml
@@ -0,0 +1,118 @@
+#%YAML
+
+#All parameters should be in SI units
+
+#Few parameters may be set to list of float or integer allowing to run successive simulation
+# there are two syntaxes:
+#       [x1, ..., xn] or range(start, end, step)
+#       eg. range(0.02, 0.09, 0.02) or [0.02, 0.04, 0.06, 0.08] will give the same results
+#       the parameter will take successively the values 0.02, 0.04, 0.06 and 0.08
+
+archi:
+    #if read_architecture is true then architecture will be constructed from the file(s) given by input_dir and input_file
+    #otherwise the architecture will be generated according to the parameters
+    read_architecture: True
+
+    #Input architecture from scanned image (distance_from_base_(mm) \t lateral_root_length_(mm) \t order)
+    #folder name
+    input_dir: ../data/
+
+    #File name:
+    #may be a list of names, eg. [file1, file2, file3] wildcar may be used
+    input_file: [plant-01.txt]
+
+    seed:
+
+    #file names with length laws relative path
+    #file format: "LR_length_mm" ; "relative_distance_to_tip"
+    #laws used to generate lateral roots of the 1st order (1_order_law), and lateral roots of order above 1 (2_order_law)
+    length_file:
+        - ../data/length*order1*.csv
+        - ../data/length*order2*.csv
+
+    #length of the primary root
+    #float or list of float
+    #unit: m
+    primary_length: 0.32
+
+    #branching delay
+    #float or list of float
+    #unit: m
+    branching_delay: 1.0e-4
+
+    #branching variability
+    #float between [0 ; 1], 0.25 means 25%
+    branching_variability: 0.25
+
+    #maximum roots order
+    order_max: 4
+
+    #vertices length
+    #unit: m
+    segment_length: 1.0e-4
+
+    #part of roots without any lateral root, distance from tip
+    #float or list of float
+    #unit: m
+    nude_length: 0.03
+
+    #reference radius of the primary root
+    #float
+    #unit: m
+    ref_radius: 7.0e-5
+
+    #radius decrease factor applied when increasing order
+    #float
+    #radius lateral order n: r = order_decrease_factor^n * ref_radius
+    order_decrease_factor: 0.7
+
+hydro:
+    #radial conductivity
+    #float
+    #unit: microL/(s.MPa.m**2)
+    k0:
+        - [0.0, 0.32]
+        - [10.0, 10.0]
+
+    #axial_conductance_data
+    #   - [x1, ......, xn]
+    #   - [K1, ....., Kn]
+    #list of float
+    #unit: microL.m/(s.Mpa)
+    axial_conductance_data:
+        - [0, 0.32]
+        - [3.0e-9,2.0e-4]
+
+experimental:
+    #water flux at the root base
+    #float
+    #unit: microL/s
+    Jv:  4.52e-3
+
+    #hydric potential outside the roots (pressure chamber)
+    #float
+    #unit: MPa
+    psi_e: 0.501325
+
+    #hydric potential at the root base (e.g. atmospheric pressure for decapitated plant)
+    #float
+    #unit: MPa
+    psi_base: .101325
+
+output:
+    #distance from the base for intercepts calculation
+    # float or list of float
+    #unit: m
+    intercepts: []
+
+    #factor to explore a k0 range
+    # float or list of float
+    radfold: 1.
+
+    #like radfold but apply to axial_conductance_data
+    axfold: 1.
+
+    #number of run with the same set of parameters i.e. number of different seeds
+    #integer
+    #enable only if read_architecture is false
+    run_nb: 1
diff --git a/src/openalea/hydroroot/analysis.py b/src/openalea/hydroroot/analysis.py
index c4718f06..0478e772 100644
--- a/src/openalea/hydroroot/analysis.py
+++ b/src/openalea/hydroroot/analysis.py
@@ -23,13 +23,13 @@ def nb_roots(g, l, root=1, dl=1e-4, max_order=None):
     """
     length = {}
 
-    if 'mylength' in g.property_names():
-        length = g.property('mylength')
+    if 'dist_to_base' in g.property_names():
+        length = g.property('dist_to_base')
     else:
         for v in pre_order2(g, root):
             pid = g.parent(v)
             length[v] = length[pid] + dl if pid else dl
-        g.properties()['mylength'] = length
+        g.properties()['dist_to_base'] = length
 
     order = None
     if max_order is not None:
diff --git a/src/openalea/hydroroot/display.py b/src/openalea/hydroroot/display.py
index 56394acd..e4c45fc6 100644
--- a/src/openalea/hydroroot/display.py
+++ b/src/openalea/hydroroot/display.py
@@ -37,9 +37,9 @@ def root_visitor(g, v, turtle, prune=prune):
         :param prune: (float) - distance from base after witch the MTG is no longer read (Default value = None)
 
         """
-        mylength = {}
-        if prune and ('mylength' in g.properties()):
-            mylength = g.property('mylength')
+        dist_to_base = {}
+        if prune and ('dist_to_base' in g.properties()):
+            dist_to_base = g.property('dist_to_base')
         angles = [90,45]+[30]*5
         n = g.node(v)
         # radius = n.radius*1.e4
@@ -49,7 +49,7 @@ def root_visitor(g, v, turtle, prune=prune):
         length = n.length * factor
 
         if prune:
-            if mylength.get(v,0.)>prune:
+            if dist_to_base.get(v,0.)>prune:
                 return
 
         if g.edge_type(v) == '+':
diff --git a/src/openalea/hydroroot/flux.py b/src/openalea/hydroroot/flux.py
index fd0fb58f..a70bac35 100644
--- a/src/openalea/hydroroot/flux.py
+++ b/src/openalea/hydroroot/flux.py
@@ -482,13 +482,13 @@ def segments_at_length(g, l, root=1, dl=1e-4):
     """
     length = {}
 
-    if 'mylength' in g.property_names():
-        length = g.property('mylength')
+    if 'dist_to_base' in g.property_names():
+        length = g.property('dist_to_base')
     else:
         for v in traversal.pre_order2(g, root):
             pid = g.parent(v)
             length[v] = length[pid] + dl if pid else dl
-        g.properties()['mylength'] = length
+        g.properties()['dist_to_base'] = length
 
     vids = []
     for v in g:
diff --git a/src/openalea/hydroroot/generator/markov.py b/src/openalea/hydroroot/generator/markov.py
index 249bfdb1..7865f8d2 100644
--- a/src/openalea/hydroroot/generator/markov.py
+++ b/src/openalea/hydroroot/generator/markov.py
@@ -309,7 +309,8 @@ def generate_g(seed = None, length_data = None, branching_variability = 0.25,
 
     # F. Bauget 2022-08-12: added if-else to be able to use the function without length data, useful for usage demo
     if length_data:
-        length_max_secondary = length_data[0].LR_length_mm.max() * 1e-3  # in m
+        #TODO avoid unit management here mm to m '* 1e-3'
+        length_max_secondary = length_data[0].iloc[:,0].max() * 1e-3  # in m
 
         law_order1 = length_law(length_data[0], scale_x = primary_length / 100., scale = segment_length)
         law_order2 = length_law(length_data[1], scale_x = length_max_secondary / 100., scale = segment_length)
diff --git a/src/openalea/hydroroot/init_parameter.py b/src/openalea/hydroroot/init_parameter.py
index bda715bc..57da49d4 100644
--- a/src/openalea/hydroroot/init_parameter.py
+++ b/src/openalea/hydroroot/init_parameter.py
@@ -200,3 +200,14 @@ def parameters_to_list(self, parameter):
                 parameter = [parameter]
 
         return parameter
+
+    def metafspm_scenario(self):
+        """
+        In metafspm scenario is just a dict
+        This function simply aggregates all variables in a single dictionary
+        :return: scenario (dict)
+        """
+        scenario = {}
+        for key, value in self.__dict__.items():
+            scenario.update(value)
+        return scenario
diff --git a/src/openalea/hydroroot/law.py b/src/openalea/hydroroot/law.py
index 3401243a..1c7101ef 100644
--- a/src/openalea/hydroroot/law.py
+++ b/src/openalea/hydroroot/law.py
@@ -232,7 +232,7 @@ def reference_relative_law(x, y, size=5e-2, scale_x=1., scale_y=1e-3):
     return length.fit_law(X, means, ext=2)
 
 
-def length_law(pd, scale_x = 1 / 100., scale_y = 1., scale = 1e-4, uniform = 'expo', size = 5):
+def length_law(pd, scale_x = 1 / 100., scale_y = 1.e-3, scale = 1e-4, uniform = 'expo', size = 5):
     """Build the function giving the lateral length according to its position on the parent branch
 
     :param pd: DataFrame
@@ -252,16 +252,16 @@ def length_law(pd, scale_x = 1 / 100., scale_y = 1., scale = 1e-4, uniform = 'ex
     * 2nd col: "relative_distance_to_tip" relative distance to tip in % so between 0 and 100.
 
     """
-    x = pd.relative_distance_to_tip.tolist()
-    y = pd.LR_length_mm.tolist()
+    x = pd.iloc[:,1].tolist()
+    y = pd.iloc[:,0].tolist()
 
     # size of the windows: in %
     size *= scale_x
-    # TODO : change '1.e-3 * scale_y' to scale_y
+
     _length_law = histo_relative_law(x, y,
                                      size = size,
                                      scale_x = scale_x,
-                                     scale_y = 1.e-3 * scale_y,
+                                     scale_y = scale_y,
                                      scale = scale,
                                      plot = False,
                                      uniform = uniform)
diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index c4759853..cdba8aad 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -407,7 +407,7 @@ def hydroroot_from_data(
 def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length = 2.0e-2, df = None, segment_length = 1.0e-4,
                   length_data = None, branching_variability = 0.25, order_max = 4.0, order_decrease_factor = 0.7,
                   ref_radius = 7.0e-5, Flag_radius = False):
-    """wrapper function: build a MTG with properties that are set like : radius, vertex length, position (distance to tip) and mylength (distance to base).
+    """wrapper function: build a MTG with properties that are set like : radius, vertex length, position (distance to tip) and dist_to_base (distance to base).
     The MTG is either generated or built from a DataFrame. The unit length in hydroroot should be in m.
 
     :param primary_length: (float) - primary root length for generated mtg (Default value = 0.13)
@@ -416,7 +416,7 @@ def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length
     :param nude_length: (float) - length from tip without lateral for generated mtg (Default value = 2.0e-2)
     :param df: (DataFrame) - DataFrame with the architecture data to be reconstructed if not None (Default value = None)
     :param segment_length: (float) - vertices length in hydroroot should be in m (Default value = 1.0e-4)
-    :param length_data: (string list) - the file name with the length data laws (Default value = None)
+    :param length_data: (list of DataFrame) - DataFrame with the length data laws (Default value = None)
     :param branching_variability: (float) (Default value = 0.25)
     :param order_max: (float) (Default value = 4.0)
     :param order_decrease_factor: (float) (Default value = 0.7)
@@ -456,16 +456,16 @@ def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length
     g = radius.compute_relative_position(g)
 
     # Calculation of the distance from base of each vertex, used for cut and flow
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['mylength'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # total_length is the total length of the RSA (sum of the length of all the segments)
     total_length = g.nb_vertices(scale = 1) * segment_length
     g, surface = radius.compute_surface(g)
-    g, volume = radius.compute_volume(g)
+    # g, volume = radius.compute_volume(g)
 
     if df is not None:
         v_base = next(g.component_roots_at_scale_iter(g.root, scale = g.max_scale()))
diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
new file mode 100644
index 00000000..fdeb50ab
--- /dev/null
+++ b/src/openalea/hydroroot/model.py
@@ -0,0 +1,565 @@
+import math
+
+from dataclasses import dataclass
+
+from openalea.mtg import traversal
+from openalea.metafspm.component import Model, declare
+from openalea.metafspm.component_factory import *
+
+from openalea.hydroroot import radius, conductance, main, flux
+from openalea.hydroroot.generator import markov
+from openalea.hydroroot.law import length_law
+from openalea.hydroroot.length import fit_law
+from openalea.hydroroot.water_solute_transport import (pressure_calculation_no_non_permeating_solutes,
+                                                       init_some_MTG_properties, pressure_calculation)
+from openalea.hydroroot.analysis import intercept
+
+EPS = 1.0e-9 # epsilon value used for convergence of Newton-Raphson loop
+
+@dataclass
+class HydroRootModel(Model):
+    # Input Parameters related to architecture
+    length_data: list = declare(default=None, unit="[adim,m]", unit_comment="relative distance to the tip, meter",
+                                description="list of float or list of 2 list of float, lateral length vs relat. dist. to tip",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    primary_length: float = declare(default=0.16, unit="m", unit_comment="", description="length of the primary root",
+                                    min_value="0", max_value="", value_comment="", references="", DOI="",
+                                    variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    seed: int = declare(default=None, unit="", unit_comment="", description="seed used to generate the architecture",
+                        min_value="0", max_value="", value_comment="", references="", DOI="",
+                        variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    branching_delay: float = declare(default=2.0e-3, unit="m", unit_comment="",
+                                     description="distance between branching laterals",
+                                     min_value="0", max_value="", value_comment="", references="", DOI="",
+                                     variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    branching_variability: float = declare(default=0.25, unit="", unit_comment="",
+                                           description="variability of the branching laterals distance",
+                                           min_value="0", max_value="1", value_comment="", references="", DOI="",
+                                           variable_type="parameter", by="HydroRootModel", state_variable_type="",
+                                           edit_by="user")
+
+    order_max: int = declare(default=4, unit="", unit_comment="", description="maximum order of laterals",
+                             min_value="0", max_value="1", value_comment="", references="", DOI="",
+                             variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
+                                    min_value="0", max_value="", value_comment="", references="", DOI="",
+                                    variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    nude_length: float = declare(default=0.02, unit="m", unit_comment="",
+                                 description="part of roots without any lateral root, distance from tip",
+                                 min_value="0", max_value="", value_comment="", references="", DOI="",
+                                 variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    ref_radius: float = declare(default=7.0e-5, unit="m", unit_comment="",
+                                description="reference radius or radius of the primary root",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    order_decrease_factor: float = declare(default=0.7, unit="", unit_comment="",
+                                           description="radius decrease factor applied when increasing order",
+                                           min_value="0", max_value="1", value_comment="", references="", DOI="",
+                                           variable_type="parameter", by="HydroRootModel", state_variable_type="",
+                                           edit_by="user")
+
+    # hydraulic parameters
+    axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+                                description="(2 list of Float) axial conductivity versus dist. to tip",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]",
+                                            unit_comment="distance to the tip, k",
+                                            description="(2 list of Float) radial conductivity versus dist. to tip",
+                                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                                            variable_type="parameter", by="HydroRootModel", state_variable_type="",
+                                            edit_by="user")
+
+    # Soil1DModel input
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                           description="Homogeneous hydrostatic potential at the vertex boundary",
+                           min_value="0", max_value="", value_comment="", references="", DOI="",
+                           variable_type="input", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # HydrostaticModel parameter
+    psi_base: float = declare(default=0.101325, unit="MPa", unit_comment="",
+                              description="Hydrostatic potential at the root base",
+                              min_value="0", max_value="", value_comment="", references="", DOI="",
+                              variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    Jv: float = declare(default=None, unit="microL.s-1", unit_comment="",
+                        description="water flux at the root base, input if invert_model=False ",
+                        min_value="0", max_value="", value_comment="", references="", DOI="",
+                        variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    invert_model: bool = declare(default=True, unit="", unit_comment="",
+                                 description="when false, distribute output flux within the root ; "
+                                             "when true, compute the output flux for the given root and conditions.",
+                                 min_value="", max_value="", value_comment="", references="", DOI="",
+                                 variable_type="parameter", by="HydrostaticModel", state_variable_type="",
+                                 edit_by="user")
+
+    cut_and_flow: bool = declare(default=False, unit="", unit_comment="",
+                                 description="Use specific model to compute conductance at tips with cut & flow.",
+                                 min_value="", max_value="", value_comment="", references="", DOI="",
+                                 variable_type="parameter", by="HydrostaticModel", state_variable_type="",
+                                 edit_by="user")
+
+    # Solute parameters
+    J_s: float = declare(default=1.0e-7, unit="mol.m-2.s-1", unit_comment="", description="Active pumping rate",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    P_s: float = declare(default=1.0e-9, unit="m.s-1", unit_comment="", description="permeability coefficient",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Cse: float = declare(default=13.96, unit="mol.m-3", unit_comment="",
+                         description="concentration of permeating solutes",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Ce: float = declare(default=0.0, unit="mol.m-3", unit_comment="",
+                        description="concentration of non-permeating solutes",
+                        min_value="0", max_value="", value_comment="", references="", DOI="",
+                        variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    sigma: float = declare(default=1.0, unit="", unit_comment="", description="reflexion coefficient",
+                           min_value="0", max_value="", value_comment="", references="", DOI="",
+                           variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    temp: float = declare(default=298.0, unit="K", unit_comment="", description="sap temperature in degree K",
+                          min_value="", max_value="", value_comment="", references="", DOI="",
+                          variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    c_base: float = declare(default=None, unit="mol.m-3", unit_comment="",
+                            description="solute concentration at the root's base",
+                            min_value="", max_value="", value_comment="", references="", DOI="",
+                            variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    # Plant scale variable
+    keq: float = declare(default=None, unit="10-9 m4.MPa-1.s-1", unit_comment="",
+                         description="equivalent conductance of the RSA",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="plant_scale_state", by="HydrostaticModel", state_variable_type="",
+                         edit_by="user")
+
+    surface: float = declare(default=None, unit="m2", unit_comment="", description="total surface of the RSA",
+                             min_value="0", max_value="", value_comment="", references="", DOI="",
+                             variable_type="plant_scale_state", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    volume: float = declare(default=None, unit="m3", unit_comment="", description="total surface of the RSA",
+                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                            variable_type="plant_scale_state", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        # self.choregrapher(module_family=self.__class__.__name__)
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+    def generate_mtg(self):
+        """
+        generate a MTG according to the input parameters using self.props
+        """
+        nb_vertices = int(self.primary_length / self.length)
+        # some distance expressed in term of nb of vertices
+        nb_nude_vertices = int(self.nude_length / self.length)
+        nb_branching_delay = int(self.branching_delay / self.length)
+
+        if isinstance(self.length_data, list):
+            length_max_secondary = self.length_data[0].iloc[:, 0]
+            law_order1 = length_law(self.length_data[0], scale_x=self.primary_length / 100., scale=self.length)
+            law_order2 = length_law(self.length_data[1], scale_x=length_max_secondary / 100., scale=self.length)
+            _length_law = [law_order1, law_order2]
+
+        self.g = markov.markov_binary_tree(nb_vertices=nb_vertices,
+                                           branching_variability=self.branching_variability,
+                                           branching_delay=nb_branching_delay,
+                                           length_law=length_law,
+                                           nude_tip_length=nb_nude_vertices,
+                                           order_max=self.order_max,
+                                           seed=self.seed)
+
+        self.compute_radius()
+        self.compute_length()
+        self.compute_dist_to_tip()
+        self.compute_surface_volume()
+        self.compute_radius()
+
+    @actual
+    def compute_metrics(self):
+        if not self.g.property('radius'):
+            # TODO ref_radius etc.
+            self.compute_radius()
+        if not self.g.property('length'):
+            self.compute_length()
+
+        self.compute_dist_to_base()
+        self.compute_dist_to_tip()
+        self.compute_surface_volume()
+        self.g.property('surface')[1] = self.surface
+        self.g.property('volume')[1] = self.volume
+
+    @actual
+    @state
+    def compute_length(self):
+        """
+        compute vertices length according to self.length
+        :return:
+        """
+        self.g = radius.compute_length(self.g, self.length)
+
+    @actual
+    @state
+    def compute_dist_to_base(self):
+        """
+        Calculation of the distance from base of each vertex, used for cut and flow
+        :return:
+        """
+        g = self.g
+        _dist_to_base = {}
+        segment_length = self.length
+        for v in traversal.pre_order2(g, 1):
+            pid = g.parent(v)
+            _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+        g.properties()['dist_to_base'] = _dist_to_base
+
+    @actual
+    @state
+    def compute_dist_to_tip(self):
+        """
+        For each vertex compute the distance from the tip of the axis bearing it in absolute
+        and relative to the length of the axis bearing it
+        :return:
+        """
+        self.g = radius.compute_relative_position(self.g)
+
+    @actual
+    @state
+    def compute_radius(self):
+        """
+        compute the radius of each vertex from self.ref_radius with fixed decrease between each order.
+        :return:
+        """
+        self.g = radius.ordered_radius(self.g, self.ref_radius, self.order_decrease_factor)
+
+    @actual
+    @state
+    def compute_surface_volume(self):
+        """
+        set self.g and compute surface, volume, etc.
+        """
+        self.g, self.surface = radius.compute_surface(self.g)
+        self.g, self.volume = radius.compute_volume(self.g)
+
+    @actual
+    @axial
+    @state
+    def compute_axial_conductance(self, data = None):
+        """
+        compute axial conductance of each vertex from self.ref_axial_conductivity
+        :return:
+        """
+        if data:
+            self.axial_conductance_data = data
+
+        if self.axial_conductance_data:
+            # Compute K using axial conductance data
+            xa, ya = self.axial_conductance_data
+            axial_conductivity_law = fit_law(xa, ya)
+
+            self.g = conductance.fit_property_from_spline(self.g, axial_conductivity_law, 'position', 'K_exp')
+            self.g = conductance.compute_K(self.g)
+
+    @actual
+    @state
+    def compute_radial_conductance(self, data = None):
+        """
+        compute radial conductance of each vertex from self.ref_radial_conductivity
+        :return:
+        """
+        if data:
+            self.k0 = data
+
+        if self.k0:
+            xr, yr = self.k0
+            radial_conductivity_law = fit_law(xr, yr)
+
+            self.g = conductance.fit_property_from_spline(self.g, radial_conductivity_law, 'position', 'k0')
+            self.g = conductance.compute_k(self.g, k0='k0')
+
+    def compute_intercepts(self, dists, dl=1e-4, max_order=None):
+        return intercept(self.g, dists, dl=dl, max_order=max_order)
+
+
+    @actual
+    @state
+    def hydrostatic_solver_flux(self):
+        """
+        Compute flux according to psi_e and psi_base
+        If self.psi_e is None use self.g.property('psi_e')
+        :return:
+        """
+        _psi_e = self.psi_e
+        f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
+        f.run()
+        self.g = f.g
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        self.keq = self.g.property('Keq')[v_base]
+        self.Jv = self.g.property('J_out')[v_base]
+
+    @actual
+    @state
+    def solute_solver_flux(self):
+        """
+        Compute flux according to psi_e and psi_base
+        If self.psi_e is None use self.g.property('psi_e')
+        :return:
+        """
+
+        # init the MTG with no solute
+        _psi_e = self.psi_e
+        f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
+        f.run()
+        self.g = f.g
+
+        self.g = init_some_MTG_properties(self.g, tau=self.J_s, Cini=self.Cse*1.0e-9, Cpeg_ini=self.Ce*1.0e-9, t=1, Ps=self.P_s)
+
+        if self.Ce <= 0.0:
+            calculation = pressure_calculation_no_non_permeating_solutes
+        else:
+            calculation = pressure_calculation
+
+        # Newton-Raphson loop
+        nb_v = self.g.nb_vertices()
+        Fdx = 1.0
+        Fdx_old = 1.
+        while Fdx > EPS:
+            g, dx, data, row, col = calculation(g, Temp=self.temp, sigma=self.sigma, Ce=self.Ce*1.0e-9, Cse=self.Cse*1.0e-9,
+                                                Pe=self.psi_e, Pbase=self.psi_base, C_base=self.c_base*1.0e-9)
+            Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
+            if abs(Fdx - Fdx_old) < EPS: break
+            Fdx_old = Fdx
+
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        self.Jv = self.g.property('J_out')[v_base]
+
+@dataclass
+class HydrostaticModel(Model):
+    # Hydraulic input from HydroRootModel
+    length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
+                                    min_value="0", max_value="", value_comment="", references="", DOI="",
+                                    variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+                                description="(2 list of Float) axial conductivity versus dist. to tip",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
+                                            description="(2 list of Float) radial conductivity versus dist. to tip",
+                                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                                            variable_type="input", by="HydroRootModel", state_variable_type="",
+                                            edit_by="user")
+
+    # Soil1DModel input
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                                description="Homogeneous hydrostatic potential at the vertex boundary",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # HydrostaticModel parameter
+    psi_base: float = declare(default=0.101325, unit="MPa", unit_comment="",
+                                description="Hydrostatic potential at the root base",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    Jv: float = declare(default=None, unit="microL.s-1", unit_comment="",
+                                description="water flux at the root base, input if invert_model=False ",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    invert_model: bool = declare(default=True, unit="", unit_comment="",
+                                description="when false, distribute output flux within the root ; "
+                                            "when true, compute the output flux for the given root and conditions.",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    cut_and_flow: bool = declare(default=False, unit="", unit_comment="",
+                                description="Use specific model to compute conductance at tips with cut & flow.",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    keq: float = declare(default=None, unit="10-9 m4.MPa-1.s-1", unit_comment="", description="equivalent conductance of the RSA",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="plant_scale_state", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+
+@dataclass
+class SoluteModel(Model):
+    # Hydraulic input from HydroRootModel
+    axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+                                description="(2 list of Float) axial conductivity versus dist. to tip",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="HydroRootModel", state_variable_type="", edit_by="user")
+
+    k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
+                                            description="(2 list of Float) radial conductivity versus dist. to tip",
+                                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                                            variable_type="input", by="HydroRootModel", state_variable_type="",
+                                            edit_by="user")
+
+    # Soil1DModel input
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                                description="Homogeneous hydrostatic potential at the vertex boundary",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # SoluteModel parameter
+    psi_base: float = declare(default=0.101325, unit="MPa", unit_comment="",
+                                description="Hydrostatic potential at the root base",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Jv: float = declare(default=None, unit="microL.s-1", unit_comment="",
+                                description="water flux at the root base, input if invert_model=False ",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    cut_and_flow: bool = declare(default=False, unit="", unit_comment="",
+                                description="Use specific model to compute conductance at tips with cut & flow.",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    J_s: float = declare(default=1.0e-7, unit="mol.m-2.s-1", unit_comment="", description="Active pumping rate",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    P_s: float = declare(default=1.0e-9, unit="m.s-1", unit_comment="", description="permeability coefficient",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Cse: float = declare(default=13.96, unit="mol.m-3", unit_comment="", description="concentration of permeating solutes",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Ce: float = declare(default=0.0, unit="mol.m-3", unit_comment="", description="concentration of non-permeating solutes",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    sigma: float = declare(default=1.0, unit="", unit_comment="", description="reflexion coefficient",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    temp: float = declare(default=298.0, unit="K", unit_comment="", description="sap temperature in degree K",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    c_base: float = declare(default=None, unit="mol.m-3", unit_comment="", description="solute concentration at the root's base",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+    def flux(self):
+        """
+        Compute flux according to psi_e and psi_base
+        If self.psi_e is None use self.g.property('psi_e')
+        :return:
+        """
+
+        # init the MTG with no solute
+        _psi_e = self.psi_e
+        f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
+        f.run()
+        self.g = f.g
+
+        self.g = init_some_MTG_properties(self.g, tau=self.J_s, Cini=self.Cse*1.0e-9, Cpeg_ini=self.Ce*1.0e-9, t=1, Ps=self.P_s)
+
+        if self.Ce <= 0.0:
+            calculation = pressure_calculation_no_non_permeating_solutes
+        else:
+            calculation = pressure_calculation
+
+        # Newton-Raphson loop
+        nb_v = self.g.nb_vertices()
+        Fdx = 1.0
+        Fdx_old = 1.
+        while Fdx > EPS:
+            g, dx, data, row, col = calculation(g, Temp=self.temp, sigma=self.sigma, Ce=self.Ce*1.0e-9, Cse=self.Cse*1.0e-9,
+                                                Pe=self.psi_e, Pbase=self.psi_base, C_base=self.c_base*1.0e-9)
+            Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
+            if abs(Fdx - Fdx_old) < EPS: break
+            Fdx_old = Fdx
+
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        self.Jv = self.g.property('J_out')[v_base]
+
+@dataclass
+class Soil1DModel(Model):
+    # Soil1DModel parameters
+    soil_data: tuple = declare(default=None, unit="[m, MPa]", unit_comment="",
+                                description="tuple of 2 lists, (z,psi_e) z=depth, psi_e=water potential",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # state variable
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                                description="Soil hydrostatic potential at the vertex boundary",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="state_variable", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+    def compute_doil_psi_e(self):
+        """
+        add a soil as heterogeneous water potential along z
+        :return:
+        """
+        self.g = main.soil_1D(self.g, self.soil_data)
diff --git a/src/openalea/hydroroot/radius.py b/src/openalea/hydroroot/radius.py
index fcf66f2a..fe217926 100644
--- a/src/openalea/hydroroot/radius.py
+++ b/src/openalea/hydroroot/radius.py
@@ -223,6 +223,9 @@ def compute_relative_position(g):
     return:
         - g
 
+    .. note::
+        For instance, needed to compute axial and radial conductances that are properties depending
+        on the distance to tip
     """
     #print 'entering MTG node positionning computation'
     scale = g.max_scale()
diff --git a/src/openalea/hydroroot/solver_wrapper.py b/src/openalea/hydroroot/solver_wrapper.py
index 4434ef8b..c5cfcef3 100644
--- a/src/openalea/hydroroot/solver_wrapper.py
+++ b/src/openalea/hydroroot/solver_wrapper.py
@@ -731,8 +731,8 @@ def Jv_cnf_calculation(g, sigma, Js, Ps):
 
         # case where the primary is shorter than laterals
         max_length = primary_length
-        mylength = g_cut[0, ig].property('mylength')
-        if max(mylength.values()) > max_length: max_length = max(mylength.values())
+        dist_to_base = g_cut[0, ig].property('dist_to_base')
+        if max(dist_to_base.values()) > max_length: max_length = max(dist_to_base.values())
 
         # set conductance
         g_cut[0, ig] = set_conductances(g_cut[0, ig], axial_pr = axial_data, k0_pr = k[0], axial_lr = axial_lr,
diff --git a/test/test_rsml_yaml.py b/test/test_rsml_yaml.py
index 949f64b6..18149a28 100644
--- a/test/test_rsml_yaml.py
+++ b/test/test_rsml_yaml.py
@@ -18,7 +18,7 @@ def root_creation(g, segment_length, ref_radius, order_decrease_factor):
     Set MTG properties and perform some gemetrical calculation
 
     The vertex radius properties is set.
-    The following properties are computed: length, position, mylength, surface, volume, total length,
+    The following properties are computed: length, position, dist_to_base, surface, volume, total length,
         primary root length
 
     :param:
@@ -40,11 +40,11 @@ def root_creation(g, segment_length, ref_radius, order_decrease_factor):
     # Calculation of the distance from base of each vertex, used for cut and flow
     # Remark: this calculation is done in flux.segments_at_length; analysis.nb_roots but there is a concern with the
     # parameter dl which should be equal to vertex length but which is not pass
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['mylength'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * segment_length
@@ -82,7 +82,7 @@ def test_rsml():
     resolution *= rsml_units_to_metre[unit]  # rsml file unit to meter
     g = import_rsml_to_discrete_mtg(g_c, segment_length = segment_length, resolution = resolution)
     
-    # calculation of g properties: radius, mylength, etc.
+    # calculation of g properties: radius, dist_to_base, etc.
     g, primary_length, _length, surface = set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor)
 
     export_mtg_to_rsml(g, "data/test_rsml_io.rsml", segment_length = segment_length)
@@ -135,7 +135,7 @@ def set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor):
     Set MTG properties and perform some gemetrical calculation
 
     The vertex radius properties is set.
-    The following properties are computed: length, position, mylength, surface, volume, total length,
+    The following properties are computed: length, position, dist_to_base, surface, volume, total length,
         primary root length
 
     :param:
@@ -157,11 +157,11 @@ def set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor):
     # Calculation of the distance from base of each vertex, used for cut and flow
     # Remark: this calculation is done in flux.segments_at_length; analysis.nb_roots but there is a concern with the
     # parameter dl which should be equal to vertex length but which is not pass
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['mylength'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * segment_length