diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 099fc74e82c08b7a4ca79330bdd9ee16102ec27f..af91ef78397535d0ec51778930df64d17b693baf 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -1150,16 +1150,16 @@ class DomainPatch(df.SubDomain):
             # assumes TP, we need to initialise interface values for the nonwetting
             # phase, too for communication, see the article.
             if self.isRichards and neighbour_assumes_TP:
-                subdom_has_phases=["wetting", "nonwetting"]
+                interface_has_phases=["wetting", "nonwetting"]
             else:
-                subdom_has_phases=self.has_phases
+                interface_has_phases=self.has_phases
 
             self.interface[interf].init_interface_values(
                 interface_marker=marker,
                 interface_marker_value=marker_value,
                 function_space=space,
                 subdomain_index=self.subdomain_index,
-                has_phases=subdom_has_phases,
+                has_phases=interface_has_phases,
                 )
 
     def _calc_gravity_expressions(self):
@@ -1200,10 +1200,14 @@ class DomainPatch(df.SubDomain):
 
         """
         t = time
+        S = self.saturation
+        saturation = dict()
         if not write_iter_for_fixed_time:
             for phase in self.has_phases:
+                saturation.update({phase: df.Function(self.function_space["pressure"]["wetting"])})
                 self.pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
                 file.write(self.pressure[phase], t)
+
             capillary_pressure = df.Function(self.function_space["pressure"]["wetting"])
             if self.isRichards:
                 capillary_pressure.assign(-self.pressure["wetting"])
@@ -1214,6 +1218,19 @@ class DomainPatch(df.SubDomain):
                 capillary_pressure.assign(pc_temp)
             capillary_pressure.rename("pc_num", "pc_num")
             file.write(capillary_pressure, t)
+
+            # write out the saturation
+            for phase in self.has_phases:
+                # saturation_w.assign(Sat_w)
+                if phase is "wetting":
+                    saturation["wetting"] = df.project(S(capillary_pressure), self.function_space["pressure"]["wetting"])
+                    saturation["wetting"].rename("Sw", "Sw")
+                else:
+                    # S_nw = 1-S(exact_capillary_pressure).vector().get_local()
+                    saturation["nonwetting"] = df.project(1-S(capillary_pressure), self.function_space["pressure"]["wetting"])
+                    saturation["nonwetting"].rename("Snw", "Snw")
+
+                file.write(saturation[phase], t)
         else:
             for phase in self.has_phases:
                 i = self.iteration_number