diff --git a/example/figs/paper_figs.py b/example/figs/paper_figs.py
index 051835f..d3e90d4 100644
--- a/example/figs/paper_figs.py
+++ b/example/figs/paper_figs.py
@@ -663,7 +663,7 @@ def plot_figure_14():
         m_df = pd.merge(obs_grouped, sims)
 
         x = m_df['west'].values + m_df['east'].values
-        y = m_df['sapWest'].values + m_df['sapEast'].values
+        y = m_df['arm2'].values + m_df['arm1'].values
 
         x_index = np.isfinite(x)
         y_index = np.isfinite(y)
@@ -674,8 +674,8 @@ def plot_figure_14():
             ax.xaxis.grid(which='minor', zorder=0)
             ax.yaxis.grid(which='major', zorder=0)
 
-            ax.plot(sims.index, sims.sapEast, c='r', linewidth=1)
-            ax.plot(sims.index, sims.sapWest, c='b', linewidth=1)
+            ax.plot(sims.index, sims['arm1'], c='r', linewidth=1)
+            ax.plot(sims.index, sims['arm2'], c='b', linewidth=1)
 
             ax.plot(obs_df.index, obs_df['east'], label='East', color='red',
                     marker='o', markeredgecolor='none', alpha=0.1)
diff --git a/example/gdc_can1_grapevine/sim.py b/example/gdc_can1_grapevine/sim.py
index ae696a3..fd657fd 100644
--- a/example/gdc_can1_grapevine/sim.py
+++ b/example/gdc_can1_grapevine/sim.py
@@ -1,5 +1,4 @@
-"""This is an example on running HydroShoot on a potted grapevine with a
-simple shoot architecture.
+"""This is an example on running HydroShoot on a Geneva Double Curtain grapevine with.
 """
 
 from os import getcwd
@@ -37,4 +36,4 @@
 # Run HydroShoot
 # =============================================================================
 
-model.run(g, str(getcwd()) + '/', scene)
+model.run(g, str(getcwd()) + '/', scene, sapflow_per_cordon=True)
diff --git a/example/gdc_can2_grapevine/sim.py b/example/gdc_can2_grapevine/sim.py
index 9567d7b..64f40db 100644
--- a/example/gdc_can2_grapevine/sim.py
+++ b/example/gdc_can2_grapevine/sim.py
@@ -35,4 +35,4 @@
 # Run HydroShoot
 # =============================================================================
 
-model.run(g, str(getcwd()) + '/', scene)
+model.run(g, str(getcwd()) + '/', scene, sapflow_per_cordon=True)
diff --git a/example/gdc_can3_grapevine/sim.py b/example/gdc_can3_grapevine/sim.py
index 58d152b..2c83b52 100644
--- a/example/gdc_can3_grapevine/sim.py
+++ b/example/gdc_can3_grapevine/sim.py
@@ -35,4 +35,4 @@
 # Run HydroShoot
 # =============================================================================
 
-model.run(g, str(getcwd()) + '/', scene)
+model.run(g, str(getcwd()) + '/', scene, sapflow_per_cordon=True)
diff --git a/src/hydroshoot/model.py b/src/hydroshoot/model.py
index 932d232..e217b97 100644
--- a/src/hydroshoot/model.py
+++ b/src/hydroshoot/model.py
@@ -16,7 +16,7 @@
 from hydroshoot.params import Params
 
 
-def run(g, wd, scene=None, write_result=True, **kwargs):
+def run(g, wd, scene=None, write_result=True, sapflow_per_cordon=False, **kwargs):
     """
     Calculates leaf gas and energy exchange in addition to the hydraulic structure of an individual plant.
 
@@ -24,6 +24,8 @@ def run(g, wd, scene=None, write_result=True, **kwargs):
     - **g**: a multiscale tree graph object
     - **wd**: string, working directory
     - **scene**: PlantGl scene
+    - **sapflow_per_cordon**: to return simulated sapflow per cordon if `True`, otherwise (default `False`) to return
+        only total sapflow of the whole plant
     - **kwargs** can include:
         - **psi_soil**: [MPa] predawn soil water potential
         - **gdd_since_budbreak**: [°Cd] growing degree-day since bubreak
@@ -100,8 +102,9 @@ def run(g, wd, scene=None, write_result=True, **kwargs):
     print 'GDD since budbreak = %d °Cd' % gdd_since_budbreak
 
     # Determination of perennial structure arms (for grapevine)
-    # arm_vid = {g.node(vid).label: g.node(vid).components()[0]._vid for vid in g.VtxList(Scale=2) if
-    #            g.node(vid).label.startswith('arm')}
+    if sapflow_per_cordon:
+        arm_vid = {g.node(vid).label: g.node(vid).components()[0]._vid for vid in g.VtxList(Scale=2) if
+                   g.node(vid).label.startswith('arm')}
 
     # Soil reservoir dimensions (inter row, intra row, depth) [m]
     soil_dimensions = params.soil.soil_dimensions
@@ -280,9 +283,11 @@ def run(g, wd, scene=None, write_result=True, **kwargs):
     # Simulations
     # ==============================================================================
 
-    sapflow = []
-    # sapEast = []
-    # sapWest = []
+    if sapflow_per_cordon:
+        sapflow = {k: [] for k in arm_vid.keys()}
+    else:
+        sapflow = {'E': []}
+
     an_ls = []
     rg_ls = []
     psi_stem = {}
@@ -358,9 +363,11 @@ def run(g, wd, scene=None, write_result=True, **kwargs):
             architecture.mtg_save(g, scene, output_path)
 
         # Plot stuff..
-        sapflow.append(g.node(vid_collar).Flux)
-        # sapEast.append(g.node(arm_vid['arm1']).Flux)
-        # sapWest.append(g.node(arm_vid['arm2']).Flux)
+        if sapflow_per_cordon:
+            for cordon in sapflow.keys():
+                sapflow[cordon].append(g.node(arm_vid[cordon]).Flux)
+        else:
+            sapflow['E'].append(g.node(vid_collar).Flux)
 
         an_ls.append(g.node(vid_collar).FluxC)
 
@@ -389,9 +396,7 @@ def run(g, wd, scene=None, write_result=True, **kwargs):
 
     # Write output
     # Plant total transpiration
-    sapflow = [flow * time_conv * 1000. for flow in sapflow]
-
-    # sapEast, sapWest = [np.array(flow) * time_conv * 1000. for i, flow in enumerate((sapEast, sapWest))]
+    sapflow = {k: [flow * time_conv * 1000. for flow in cordon_flow] for k, cordon_flow in sapflow.items()}
 
     # Median leaf temperature
     t_ls = [np.median(Tlc_dict[date].values()) for date in meteo.time]
@@ -402,12 +407,11 @@ def run(g, wd, scene=None, write_result=True, **kwargs):
     results_dict = {
         'Rg': rg_ls,
         'An': an_ls,
-        'E': sapflow,
-        # 'sapEast': sapEast,
-        # 'sapWest': sapWest,
         'Tleaf': t_ls
     }
 
+    results_dict.update(sapflow)
+
     # Results DataFrame
     results_df = DataFrame(results_dict, index=meteo.time)