diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index af91ef78397535d0ec51778930df64d17b693baf..217e4ece1650773482fa90fa9a149f9328ae5b55 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -504,6 +504,7 @@ class DomainPatch(df.SubDomain):
 
                 flux_function = df.project(flux, self.function_space["flux"][phase])
                 flux_x, flux_y = flux_function.split()
+                # functions get initialised with values zero.
                 gl0 = df.Function(self.function_space["gli"][phase])
                 # # get dictionaries of global dof indices and coordinates along
                 # # the interface. These dictionaries have the facet indices of
@@ -528,7 +529,9 @@ class DomainPatch(df.SubDomain):
                         flux_x_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_x"]
                         flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
 
-                        if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
+                        if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+                            # note that even if we include gravity, the correct
+                            # flux has been set above (flux=0 without gravity).
                             gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
                                 + flux_y.vector()[flux_y_dof_index]*normal[1]
                         else: