diff --git a/srcSubcycle/AMRLevelMushyLayer.H b/srcSubcycle/AMRLevelMushyLayer.H
index bc1f86f..55c7b49 100644
--- a/srcSubcycle/AMRLevelMushyLayer.H
+++ b/srcSubcycle/AMRLevelMushyLayer.H
@@ -275,6 +275,9 @@ public:
   /// Initial data for a pure fluid with top and bottom cold boundaries
   void initialDataLens();
 
+  /// Initial data for a linear enthalpy gradient
+  void initialDataLinearGradient();
+
   /// Initial data for a test case with regions of zero porosity
   void initialDataZeroPorosityTest();
 
diff --git a/srcSubcycle/AMRLevelMushyLayer.cpp b/srcSubcycle/AMRLevelMushyLayer.cpp
index 6b09efa..ee2fa83 100644
--- a/srcSubcycle/AMRLevelMushyLayer.cpp
+++ b/srcSubcycle/AMRLevelMushyLayer.cpp
@@ -786,46 +786,47 @@ Real AMRLevelMushyLayer::advance()
   return -1;
 }
 
-void AMRLevelMushyLayer::addHeatSource(LevelData<FArrayBox> &src)
+void AMRLevelMushyLayer::addHeatSource(LevelData<FArrayBox>& src)
 {
   // This is where we add in a heat source, if required
-  Real gaussian_heat_source_size = 0.0;
-  Real gaussian_heat_source_width = 0.0;
-  Real gaussian_heat_source_depth = 0.0;
-  Real gaussian_heat_source_xpos = m_domainWidth / 2;
-
-  ParmParse ppHeatSource("heatSource");
-  ppHeatSource.query("size", gaussian_heat_source_size);
-  ppHeatSource.query("width", gaussian_heat_source_width);
-  ppHeatSource.query("depth", gaussian_heat_source_depth);
+      Real gaussian_heat_source_size = 0.0;
+      Real gaussian_heat_source_width = 0.0;
+      Real gaussian_heat_source_depth = 0.0;
+      Real gaussian_heat_source_xpos = m_domainWidth/2;
+
+      ParmParse ppHeatSource("heatSource");
+      ppHeatSource.query("size", gaussian_heat_source_size);
+      ppHeatSource.query("width", gaussian_heat_source_width);
+      ppHeatSource.query("depth", gaussian_heat_source_depth);
   //    ppHeatSource.query("depth", gaussian_heat_source_depth);
-  ppHeatSource.query("xpos", gaussian_heat_source_xpos);
+      ppHeatSource.query("xpos", gaussian_heat_source_xpos);
 
-  if (gaussian_heat_source_size != 0.0)
-  {
-    for (DataIterator dit = m_grids.dataIterator(); dit.ok(); ++dit)
-    {
-      for (BoxIterator bit = BoxIterator(m_grids[dit]); bit.ok(); ++bit)
+      // enthalpy is in the first component of the source term
+      int Hcomp = 0;
+
+      if (gaussian_heat_source_size != 0.0)
       {
-        IntVect iv = bit();
-        RealVect loc;
-        ::getLocation(iv, loc, m_dx);
+        for (DataIterator dit = m_grids.dataIterator(); dit.ok(); ++dit)
+        {
+          for (BoxIterator bit = BoxIterator(m_grids[dit]); bit.ok(); ++bit)
+          {
+            IntVect iv = bit();
+            RealVect loc;
+            ::getLocation(iv, loc, m_dx);
 
-        // Set the 0th component (enthalpy)
-        src[dit](iv, 0) =
-            gaussian_heat_source_size /
-            (gaussian_heat_source_width * sqrt(2 * M_PI)) *
-            exp(-0.5 * pow((loc[0] - gaussian_heat_source_xpos) /
-                               gaussian_heat_source_width,
-                           2)) *
-            0.5 *
-            (1 + tanh(10 * (loc[1] -
-                            (m_domainHeight - gaussian_heat_source_depth))));
+            Real porosity = (*m_scalarNew[ScalarVars::m_porosity])[dit](iv);
+
+            src[dit](iv, Hcomp) = (1-porosity)*gaussian_heat_source_size/(gaussian_heat_source_width*sqrt(2*M_PI))
+                * exp(-0.5*pow((loc[0]-gaussian_heat_source_xpos)/gaussian_heat_source_width, 2))
+                * 0.5*(1 + tanh(10*(loc[1]-(m_domainHeight-gaussian_heat_source_depth) ) ));
+
+          }
+
+        }
       }
-    }
-  }
 }
 
+
 void AMRLevelMushyLayer::setVelZero(FArrayBox &a_vel,
                                     const FArrayBox &a_porosity,
                                     const Real a_limit, const int a_radius)
diff --git a/srcSubcycle/AMRLevelMushyLayerInit.cpp b/srcSubcycle/AMRLevelMushyLayerInit.cpp
index 3110538..6633648 100644
--- a/srcSubcycle/AMRLevelMushyLayerInit.cpp
+++ b/srcSubcycle/AMRLevelMushyLayerInit.cpp
@@ -1230,6 +1230,30 @@ void AMRLevelMushyLayer::initialDataLens()
 
 }
 
+void AMRLevelMushyLayer::initialDataLinearGradient()
+{
+  // Iterate over the different boxes on this level of refinement
+  DataIterator dit = m_grids.dataIterator();
+  for (dit.reset(); dit.ok(); ++dit)
+  {
+    // Iterate over the cells within each box
+    BoxIterator bit((*m_scalarNew[0])[dit].box());
+    for (bit.reset(); bit.ok(); ++bit)
+    {
+      // This is the location of the grid cell in indexing space
+      IntVect iv = bit();
+      // Get the location of this grid cell in x-y space
+      RealVect loc;
+      getLocation(iv, loc, m_dx);
+      // Set the enthalpy = -stefan number*z, where z is the vertical co-ordinate
+      // note that loc[1] is the vertical co-ordinate
+      (*m_scalarNew[ScalarVars::m_enthalpy])[dit](iv) = -m_parameters.stefan*loc[1];
+      // Set the bulk concentration
+      (*m_scalarNew[ScalarVars::m_bulkConcentration])[dit](iv) = -1.0;
+    }
+  }
+}
+
 void AMRLevelMushyLayer::initialDataSidewallHeating()
 {
   // heating in x direction, i.e. left/right walls;
@@ -2139,6 +2163,9 @@ void AMRLevelMushyLayer::initialData()
     case 2:
       initialDataLens();
       break;
+    case 3:   // if main.initData = 2, use this function:
+      initialDataLinearGradient();
+      break;
     default:
       initialDataDefault();
     }