diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp
index 26f734efd64..f5e918d2e02 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp
@@ -139,6 +139,14 @@ void SourceFluxStatsAggregator::registerDataOnMesh( Group & meshBodies )
       } );
     }
   } );
+
+  if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 )
+  {
+    std::ofstream outputFile( m_csvFilename );
+    TableCSVFormatter const tableStatFormatter( m_csvLayout );
+    outputFile << tableStatFormatter.headerToString();
+    outputFile.close();
+  }
 }
 
 void SourceFluxStatsAggregator::gatherStatsForLog( bool logLevelActive,
@@ -215,15 +223,15 @@ void SourceFluxStatsAggregator::outputStatsToCSV( TableData & csvData )
 {
   if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 )
   {
-    std::ofstream outputFile( m_csvFilename );
+    std::ofstream outputFile( m_csvFilename, std::ios::app );
     TableCSVFormatter const tableStatFormatter( m_csvLayout );
-    outputFile << tableStatFormatter.toString( csvData );
+    outputFile << tableStatFormatter.dataToString( csvData );
     outputFile.close();
     csvData.clear();
   }
 }
 
-bool SourceFluxStatsAggregator::execute( real64 const GEOS_UNUSED_PARAM( time_n ),
+bool SourceFluxStatsAggregator::execute( real64 const GEOS_UNUSED_PARAM ( time_n ),
                                          real64 const GEOS_UNUSED_PARAM( dt ),
                                          integer const GEOS_UNUSED_PARAM( cycleNumber ),
                                          integer const GEOS_UNUSED_PARAM( eventCounter ),
@@ -241,30 +249,31 @@ bool SourceFluxStatsAggregator::execute( real64 const GEOS_UNUSED_PARAM( time_n
   {
     TableData logData;
     TableData csvData;
-    meshLevelStats.stats() = StatData();
+    meshLevelStats.stats().reset();
     forAllFluxStatsWrappers( meshLevel,
                              [&] ( MeshLevel &, WrappedStats & fluxStats )
     {
-      fluxStats.stats() = StatData();
+      fluxStats.stats().reset();
       forAllRegionStatsWrappers( meshLevel, fluxStats.getFluxName(),
                                  [&] ( ElementRegionBase & region, WrappedStats & regionStats )
       {
-        regionStats.stats() = StatData();
+        regionStats.stats().reset();
 
         forAllSubRegionStatsWrappers( region, regionStats.getFluxName(),
                                       [&] ( ElementSubRegionBase &, WrappedStats & subRegionStats )
         {
+          subRegionStats.stats().reset();
           subRegionStats.finalizePeriod();
-          regionStats.stats().combine( subRegionStats.stats() );
+          regionStats.combine( subRegionStats );
         } );
-        fluxStats.stats().combine( regionStats.stats() );
+        fluxStats.combine( regionStats );
 
         gatherStatsForLog( regionsStatsOn,
                            fluxStats.getFluxName(), region.getName(), logData, regionStats );
         gatherStatsForCSV( fluxStats.getFluxName(), region.getName(), csvData, regionStats );
       } );
 
-      meshLevelStats.stats().combine( fluxStats.stats() );
+      meshLevelStats.combine( fluxStats );
 
       gatherStatsForLog( fluxesStatsOn,
                          fluxStats.getFluxName(), allRegionsStr, logData, fluxStats );
@@ -370,7 +379,7 @@ void SourceFluxStatsAggregator::WrappedStats::finalizePeriod()
 
   // produce the period stats of this rank
   m_stats.m_elementCount = m_periodStats.m_elementCount;
-  m_statsPeriodStart = m_periodStats.m_periodStart;
+  m_statsPeriodStart = MpiWrapper::max( m_periodStats.m_periodStart );
   m_statsPeriodDT = m_periodStats.m_timeStepDeltaTime + m_periodStats.m_periodPendingDeltaTime;
 
   real64 const timeDivisor = m_statsPeriodDT > 0.0 ? 1.0 / m_statsPeriodDT : 0.0;
@@ -387,6 +396,11 @@ void SourceFluxStatsAggregator::WrappedStats::finalizePeriod()
   // start a new timestep
   m_periodStats.reset();
 }
+void SourceFluxStatsAggregator::WrappedStats::combine( WrappedStats const & other )
+{
+  stats().combine( other.stats() );
+  m_statsPeriodStart = LvArray::math::max( m_statsPeriodStart, other.m_statsPeriodStart );
+}
 void SourceFluxStatsAggregator::WrappedStats::PeriodStats::allocate( integer phaseCount )
 {
   if( m_timeStepMass.size() < phaseCount )
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp
index b777cd3ba07..5b015e22232 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp
@@ -112,6 +112,14 @@ class SourceFluxStatsAggregator final : public FieldStatisticsBase< FlowSolverBa
      */
     void finalizePeriod();
 
+    /**
+     * @brief Aggregate the statistics of the instance with those of another one.
+     * @details Combines the statistics from the other object into this one and also advances the period
+     * start if the other object has a later time recorded.
+     * @param other the other WrappedStats object.
+     */
+    void combine( WrappedStats const & other );
+
     /**
      * @return the reference to the wrapped stats data collected over the last period (one timestep or more), computed by finalizePeriod()
      */
@@ -151,7 +159,7 @@ class SourceFluxStatsAggregator final : public FieldStatisticsBase< FlowSolverBa
     /// stats data collected over the last period (one timestep or more), computed by finalizePeriod()
     StatData m_stats;
     /// the start time of the wrapped stats period (in s)
-    real64 m_statsPeriodStart;
+    real64 m_statsPeriodStart{-LvArray::NumericLimits< real64 >::max};
     /// the duration of the wrapped stats period (in s)
     real64 m_statsPeriodDT;
 
@@ -169,7 +177,7 @@ class SourceFluxStatsAggregator final : public FieldStatisticsBase< FlowSolverBa
       /// time that the current timestep is simulating.
       real64 m_timeStepDeltaTime = 0.0;
       /// start time of the current period.
-      real64 m_periodStart = 0.0;
+      real64 m_periodStart = -LvArray::NumericLimits< real64 >::max;
       /// delta time from all previous time-step of the current period.
       real64 m_periodPendingDeltaTime = 0.0;
       /// number of cell elements targeted by this instance