diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 8cb8c737e0993e7ca25f984327bd4b83c08734f2..146ee4634e5730dfb7d4bd06d5c2c530e88437df 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -281,7 +281,7 @@ class LDDsimulation(object):
         self.gravity_acceleration = gravity_acceleration
         self.densities = densities
         if self.include_gravity and self.densities is None:
-                raise RuntimeError("include_gravity = True but no densities were given.")
+            raise RuntimeError("include_gravity = True but no densities were given.")
         self._parameters_set = True
 
 
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 3c4572a9a167f760af68cd4f5fd6dc7b18a75511..240db36661a0a3744421c0a4faf094962c865b82 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -392,27 +392,23 @@ class DomainPatch(df.SubDomain):
             form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
             rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
             if self.include_gravity:
-                # z_a included the minus sign already, see self._calc_gravity_expressions
-                rhs_gravity = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
+                rhs_gravity = (-dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a)))*dx
         else:
             form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
             # the non wetting phase has a + before instead of a - before the bracket
             rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
             if self.include_gravity:
-                # z_a included the minus sign, see self._calc_gravity_expressions
-                rhs_gravity = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
+                rhs_gravity = (-dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a)))*dx
 
         if self.has_interface is not None:
             form = form1 + form2 + sum(interface_forms)
             if self.include_gravity:
-                # z_a included the minus sign already, see self._calc_gravity_expressions
                 rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
             else:
                 rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
         else:
             form = form1 + form2
             if self.include_gravity:
-                # z_a included the minus sign already, see self._calc_gravity_expressions
                 rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx + rhs_gravity
             else:
                 rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx
@@ -478,8 +474,7 @@ class DomainPatch(df.SubDomain):
                     # the wetting phase is always present and we always need
                     # to calculate a gl0 term.
                     if self.include_gravity:
-                        # z_a included the minus sign, see self._calc_gravity_expressions
-                        flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev - z_a)
+                        flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
                     else:
                         flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
                 else:
@@ -487,12 +482,12 @@ class DomainPatch(df.SubDomain):
                     if self.isRichards: # and phase is "nonwetting"
                         # in this case the current interface has a TP neighbour
                         if self.include_gravity:
-                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
+                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(z_a)
                         else:
                             flux = 0
                     else:
                         if self.include_gravity:
-                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
+                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
                         else:
                             flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
 
@@ -1171,7 +1166,7 @@ class DomainPatch(df.SubDomain):
         for phase in has_phases:
             rho = df.Constant(self.densities[phase])
             self.gravity_term.update(
-                {phase: df.Expression('-g*rho*x[1]',
+                {phase: df.Expression('g*rho*x[1]',
                                       domain=self.mesh, #
                                       degree=self.interpolation_degree,
                                       g=g,