diff --git a/bin/HEAT1D/graphFunctions.cpp b/bin/HEAT1D/graphFunctions.cpp
index ebb549f0..634db2ce 100644
--- a/bin/HEAT1D/graphFunctions.cpp
+++ b/bin/HEAT1D/graphFunctions.cpp
@@ -18,7 +18,8 @@ outputType getOutputType(outputGroup myOut)
 
         return outputType::profile;
     else if (myOut == energyBalance ||
-             myOut == errorBalance)
+             myOut == errorBalance ||
+             myOut == bottom)
         return outputType::single;
     else
         return outputType::single;
@@ -50,6 +51,8 @@ void setColorScale(Crit3DColorScale* myColorScale, outputGroup outGroup, Crit3DO
             myColorScale->setMaximum(3);
         else if (outGroup == outputGroup::surfaceResistances)
             myColorScale->setMaximum(2);
+        else if (outGroup == outputGroup::bottom)
+            myColorScale->setMaximum(1);
 
         *graphLinear = true;
     }
@@ -100,6 +103,9 @@ QVector<QPointF> getSingleSeries(Crit3DOut* myOut, outputVar myVar, float* MINVA
             case outputVar::MBR_water :
                 myVal = myOut->errorOutput[i].waterMBR.y();
                 break;
+
+            case outputVar::bottomDrainage :
+                myVal = myOut->bottomFluxes[i].drainage.y();
         }
 
         myPoint.setY(myVal);
diff --git a/bin/HEAT1D/graphFunctions.h b/bin/HEAT1D/graphFunctions.h
index 2f5526b3..cab8a825 100644
--- a/bin/HEAT1D/graphFunctions.h
+++ b/bin/HEAT1D/graphFunctions.h
@@ -23,7 +23,8 @@ enum outputGroup
     energyBalance = 10,
     surfaceResistances = 11,
     soilHeatConductivity = 12,
-    errorBalance = 13
+    bottom = 13,
+    errorBalance = 14
 };
 
 enum outputVar
@@ -34,7 +35,8 @@ enum outputVar
     aerodynamicResistence = 3,
     soilSurfaceResistence = 4,
     MBR_heat = 5,
-    MBR_water = 6
+    MBR_water = 6,
+    bottomDrainage = 7
 };
 
 static const char * outputVarString[] =
@@ -84,6 +86,10 @@ struct heatErrors{
     QPointF waterMBR;
 };
 
+struct bottomFlux{
+    QPointF drainage;
+};
+
 class Crit3DOut
 {
 public:
@@ -100,6 +106,7 @@ class Crit3DOut
     QVector<landSurfaceStatus> landSurfaceOutput;
     QVector<profileStatus> profileOutput;
     QVector<heatErrors> errorOutput;
+    QVector<bottomFlux> bottomFluxes;
 
     QString getTextOutput(outputGroup myOutGroup);
 };
diff --git a/bin/HEAT1D/heat1D.cpp b/bin/HEAT1D/heat1D.cpp
index 3e5e3a31..2c0b2bb5 100644
--- a/bin/HEAT1D/heat1D.cpp
+++ b/bin/HEAT1D/heat1D.cpp
@@ -7,9 +7,6 @@
 
 #include "heat1D.h"
 #include "soilFluxes3D.h"
-#include "types.h"
-#include "boundary.h"
-#include "soilPhysics.h"
 #include "commonConstants.h"
 
 Crit3DOut output;
@@ -261,6 +258,7 @@ bool initializeHeat1D(bool useInputSoils)
 
             result = soilFluxes3D::setWaterContent(indexNode, ((indexNode-1)*(myHeat1D.initialSaturationBottom - myHeat1D.initialSaturationTop)
                                                                    / (myHeat1D.NodesNumber-2)+myHeat1D.initialSaturationTop) * (myThetaS - myThetaR) + myThetaR);
+
             if (result != CRIT3D_OK) printf("\n error in SetWaterContent!");
 
             if (myHeat1D.computeHeat)
@@ -336,11 +334,13 @@ void getOutputAllPeriod(long firstIndex, long lastIndex, Crit3DOut *output)
     profileStatus myProfile;
     landSurfaceStatus mySurfaceOutput;
     heatErrors myErrors;
+    bottomFlux myBottomFluxes;
 
     output->nrValues++;
     output->profileOutput.push_back(myProfile);
     output->landSurfaceOutput.push_back(mySurfaceOutput);
     output->errorOutput.push_back(myErrors);
+    output->bottomFluxes.push_back(myBottomFluxes);
 
     for (myIndex = firstIndex ; myIndex <= lastIndex ; myIndex++ )
     {
@@ -460,6 +460,11 @@ void getOutputAllPeriod(long firstIndex, long lastIndex, Crit3DOut *output)
     myValue = soilFluxes3D::getWaterMBR();
     if (isValid(myValue)) myPoint.setY(myValue);
     output->errorOutput[output->nrValues-1].waterMBR = myPoint;
+
+    //bottom fluxes
+    myValue = soilFluxes3D::getBoundaryWaterFlow(lastIndex);
+    if (isValid(myValue)) myPoint.setY(myValue);
+    output->bottomFluxes[output->nrValues-1].drainage = myPoint;
 }
 
 QString Crit3DOut::getTextOutput(outputGroup outGroup)
diff --git a/bin/HEAT1D/linearGraph.cpp b/bin/HEAT1D/linearGraph.cpp
index eca05a1d..c5aace2e 100644
--- a/bin/HEAT1D/linearGraph.cpp
+++ b/bin/HEAT1D/linearGraph.cpp
@@ -143,6 +143,14 @@ void LinearGraph::setProperties(outputGroup graphType)
         startIndex = 5;
 
     }
+    else if (graphType == outputGroup::bottom)
+    {
+        axisY->setTitleText(QString("water flow [%1]").arg(QString("m3")));
+
+        myCurveNames.append("water flow");
+        startIndex = 7;
+
+    }
 }
 
 void LinearGraph::tooltipLineSeries(QPointF point, bool state)
diff --git a/bin/HEAT1D/mainwindow.cpp b/bin/HEAT1D/mainwindow.cpp
index 68644974..de31d32b 100644
--- a/bin/HEAT1D/mainwindow.cpp
+++ b/bin/HEAT1D/mainwindow.cpp
@@ -52,6 +52,7 @@ MainWindow::MainWindow(QWidget *parent) :
     ui->listWidget->addItem("surface energy balance (W m-2)");
     ui->listWidget->addItem("surface resistances (s m-1)");
     ui->listWidget->addItem("heat conductivity (W m-1 K-1)");
+    ui->listWidget->addItem("bottom fluxes (m3)");
     ui->listWidget->addItem("error balance ()");
 }
 
diff --git a/bin/HEAT1D/mainwindow.ui b/bin/HEAT1D/mainwindow.ui
index 5bf33a3a..6426d725 100644
--- a/bin/HEAT1D/mainwindow.ui
+++ b/bin/HEAT1D/mainwindow.ui
@@ -51,7 +51,7 @@
      </rect>
     </property>
     <property name="currentIndex">
-     <number>1</number>
+     <number>0</number>
     </property>
     <widget class="QWidget" name="tab">
      <attribute name="title">