diff --git a/domainPatch.py b/domainPatch.py
index 3855a6d960a0abdbb144302a884edb76c5c1fcea..c92816d5d527b4759b1f74c9f2baca076986d4dc 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -196,13 +196,17 @@ class DomainPatch(df.SubDomain):
         ds = self.ds
         dt = self.timestep_size
         if self.isRichards:
+            if phase == 'nonwetting':
+                print("CAUTION: You invoked domainPatch.governing_problem() with \n#
+                       Parameter phase = 'nonwetting'. But isRichards = True for \n#
+                       current subdomain. \nReturning wetting phase equation only.\n")
             Lw = self.L['wetting']
             # this is pw_i in the paper
             pw = self.pressure['wetting']
             # this is pw_{i-1} in the paper
             pw_prev_iter = self.pressure_prev_iter['wetting']
             pw_prev_timestep = self.pressure_prev_timestep['wetting']
-            v  = self.testfunction['wetting']
+            phi_w  = self.testfunction['wetting']
             Lambda = self.lambda_param['wetting']
             kw = self.relative_permeability['wetting']
             S = self.saturation
@@ -210,30 +214,101 @@ class DomainPatch(df.SubDomain):
             source_w = self.source['wetting']
             # we need to have all interfaces in the form
             interface_forms = []
-            for index in self.has_interface:
-                interface_forms.append((dt*Lambda*pw*v)*ds[index])
-            form = (Lw*pw*v)*dx + (dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v)))*dx + df.sum(interface_forms)
-            # assemble rhs for
-            rhs_interface_forms = []
-            self._calc_gli_term(iteration = iter_num)
-            # gli = self.gli['wetting']
-            for index in self.has_interface:
-                gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
-                rhs_interface_forms.append((dt*gli*v)*ds[index])
-            rhs = (Lw*pw_prev_iter*v)*dx - ((S(pw_prev_iter) - S(pw_prev_timestep)*v)*dx + dt*source_w*v)*dx - df.sum(rhs_interface_forms)
+            for interface in self.has_interface:
+                interface_forms.append((dt*Lambda*pw*phi_w)*ds[interface])
+            form1 = (Lw*pw*phi_w)*dx
+            form2 = dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v))*dx
+            form = form1 + form2 + df.sum(interface_forms)
+            # # assemble rhs for
+            # rhs_interface_forms = []
+            # self._calc_gli_term(iteration = iter_num)
+            # # gli = self.gli['wetting']
+            # for index in self.has_interface:
+            #     gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
+            #     rhs_interface_forms.append((dt*gli*phi_w)*ds[index])
+            rhs1 = (Lw*pw_prev_iter*phi_w)*dx
+            rhs2 = -((S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
+            rhs_without_gli = rhs1 + rhs2 + dt*source_w*phi_w*dx
             form_and_rhs = {#
                 'form': form,#
-                'rhs' : rhs
+                'rhs_without_gli' : rhs_without_gli
             }
+            return form_and_rhs
         else:
+            La = self.L['alpha']
+            # Lnw = self.L['nonwetting']
+            # this is pw_i/pnw_i in the paper
+            pa = self.pressure['alpha']
+            # pnw = self.pressure['nonwetting']
+            # this is pw_{i-1} in the paper
+            pw_prev_iter = self.pressure_prev_iter['wetting']
+            pnw_prev_iter = self.pressure_prev_iter['nonwetting']
+            pw_prev_timestep = self.pressure_prev_timestep['wetting']
+            pnw_prev_timestep = self.pressure_prev_timestep['nonwetting']
+            pc_prev_iter = pnw_prev_iter - pw_prev_iter
+            pc_prev_timestep = pnw_prev_timestep - pw_prev_timestep
+            phi_a  = self.testfunction['alpha']
+            # phi_nw  = self.testfunction['nonwetting']
+            Lambda_a = self.lambda_param['alpha']
+            # Lambda_nw = self.lambda_param['nonwetting']
+            ka = self.relative_permeability['alpha']
+            # knw = self.relative_permeability['nonwetting']
+            S = self.saturation
+            mu_a = self.viscosity['alpha']
+            # mu_nw = self.viscosity['nonwetting']
+            source_a = self.source['alpha']
+            # source_nw = self.source['nonwetting']
+
             if phase == "wetting":
-                print("governing form wetting phase")
+                # gather interface forms
+                interface_forms = []
+                for interface in self.has_interface:
+                    interface_forms.append((dt*Lambda_a*pa*phi_a)*ds[interface])
+
+                form1 = (La*pa*phi_a)*dx
+                form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
+                form = form1 + form2 + df.sum(interface_forms)
+
+                rhs1 = (La*pw_prev_iter*phi_a)*dx
+                rhs2 = -((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
+                rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a)*dx
+
+                form_and_rhs = {#
+                    'form': form,#
+                    'rhs_without_gli' : rhs_without_gli
+                }
+                return form_and_rhs
+
             elif phase == "nonwetting":
-                print("governing form nonwetting phase")
-            else:
-                raise RuntimeError('missmatch for input parameter phase. \nExpected either phase == wetting or phase == nonwetting ')
+                # gather interface forms
+                interface_forms = []
+                for interface in self.has_interface:
+                    if self.interface[interface].neighbour_is_Richards(self.subdomain_index):
+                        # if the neighbour of our subdomain (with index self.subdomain_index)
+                        # assumes a Richards model, we assume constant normalised atmospheric pressure
+                        # and no interface term is practically appearing.
+                        interface_forms.append((df.Constant(0)*phi_a)*ds[interface])
+                    else:
+                        interface_forms.append((dt*Lambda_a*pa*phi_a)*ds[interface])
+
+                form1 = (La*pa*phi_a)*dx
+                form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
+                form = form1 + form2 + df.sum(interface_forms)
+
+                rhs1 = (La*pnw_prev_iter*phi_a)*dx
+                # the non wetting phase has a + before instead of a - before the bracket
+                rhs2 = ((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
+                rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a)*dx
+
+                form_and_rhs = {#
+                    'form': form,#
+                    'rhs_without_gli' : rhs_without_gli
+                }
+                return form_and_rhs
 
-        return form
+            else:
+                raise RuntimeError('missmatch for input parameter phase. \n#
+                Expected either phase == wetting or phase == nonwetting ')
 
     #### PRIVATE METHODS
     def _init_function_space(self) -> None:
@@ -357,103 +432,97 @@ class DomainPatch(df.SubDomain):
                      'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
                 )
 
-    # def _calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
-    #     """calculate the gl terms for the iteration number iteration
-    #
-    #     calculate the gl terms for the iteration number iteration and save them
-    #     to self.gli.
-    #     """
-    #     for ind in self.has_interface:
-    #         # self._calc_gli_term_on(interface, iteration)
-    #         interface = self.interface[ind]
-    #         subdomain = self.subdomain_index
-    #         # neighbour index
-    #         neighbour = interface.neighbour[subdomain]
-    #         neighbour_iter_num = interface.current_iteration[neighbour]
-    #
-    #         Lambda = self.lambda_param
-    #         if iteration == 0:
-    #             n = df.FacetNormal(self.mesh)
-    #             pw_prev_iter = self.pressure_prev_timestep
-    #             k = self.relative_permeability
-    #             S = self.saturation
-    #             mu = self.viscosity
-    #
-    #             if not neighbour_iter_num == 0:
-    #                 # save the gl term from the neighbour first to use it in the
-    #                 # next iteration
-    #                 if self.isRichards:
-    #                     interface.gli_term_prev[subdomain].update(#
-    #                         {'wetting': interface.gli_term[neighbour]['wetting']}
-    #                         )
-    #                 else:
-    #                     interface.gli_term_prev[subdomain].update(#
-    #                         {'wetting': interface.gli_term[neighbour]['wetting']},#
-    #                         {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
-    #                         )
-    #
-    #             if self.isRichards:
-    #                 mu_w = mu['wetting']
-    #                 kw = k['wetting']
-    #                 pwn_prev = pw_prev_iter['wetting']
-    #                 #calc gl0 from the flux term
-    #                 flux = -1/mu_w*df.dot(kw(S(pwn_prev))*df.grad(pwn_prev)
-    #                 gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
-    #                 interface.gli_term[subdomain].update({'wetting': gl0})
-    #             else:
-    #                 mu_w = mu['wetting']
-    #                 mu_nw = mu['nonwetting']
-    #                 pwn_prev = pw_prev_iter['wetting']
-    #                 pnwn_prev = pw_prev_iter['nonwetting']
-    #                 kw = k['wetting']
-    #                 knw = k['nonwetting']
-    #                 #calc gl0 from the flux term
-    #                 fluxw = -1/mu_w*df.dot(kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
-    #                 fluxnw = -1/mu_nw*df.dot(knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
-    #                 gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
-    #                 gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
-    #                 interface.gli_term[subdomain].update({'wetting': gwl0})
-    #                 interface.gli_term[subdomain].update({'nonwetting': gnwl0})
-    #
-    #             interface.current_iteration[subdomain] += 1
-    #
-    #         else:
-    #             if neighbour_iter_num == iteration:
-    #                 if self.isRichards:
-    #                     interface.gli_term_prev[subdomain].update(#
-    #                         {'wetting': interface.gli_term[neighbour]['wetting']}
-    #                         )
-    #
-    #                     pwn_prev = pw_prev_iter['wetting']
-    #                     #calc gl0 from the flux term
-    #                     flux = -1/mu_w*df.dot(kw(S(pwn_prev))*df.grad(pwn_prev)
-    #                     gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
-    #                     interface.gli_term[subdomain].update({'wetting': gl0})
-    #                 else:
-    #                     mu_w = mu['wetting']
-    #                     mu_nw = mu['nonwetting']
-    #                     pwn_prev = pw_prev_iter['wetting']
-    #                     pnwn_prev = pw_prev_iter['nonwetting']
-    #                     kw = k['wetting']
-    #                     knw = k['nonwetting']
-    #                     #calc gl0 from the flux term
-    #                     fluxw = -1/mu_w*df.dot(kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
-    #                     fluxnw = -1/mu_nw*df.dot(knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
-    #                     gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
-    #                     gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
-    #                     interface.gli_term[subdomain].update({'wetting': gwl0})
-    #                     interface.gli_term[subdomain].update({'nonwetting': gnwl0})
-    #
-    #             else:
-    #                 print('')
-    #             interface.current_iteration[subdomain] += 1
-    # # def _calc_gli_term_on(self, interface: int, iteration: int) -> None:
-    # #     """calculate the gl terms for the iteration number iteration
-    # #
-    # #     calculate the gl terms for the iteration number iteration and save them
-    # #     to self.gli.
-    # #     """
-    # # END is_Richards
+    def _calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
+        """calculate the gl terms for the iteration number iteration
+
+        calculate the gl terms for the iteration number iteration and save them
+        to self.gli.
+        """
+        for ind in self.has_interface:
+            # self._calc_gli_term_on(interface, iteration)
+            interface = self.interface[ind]
+            subdomain = self.subdomain_index
+            # neighbour index
+            neighbour = interface.neighbour[subdomain]
+            neighbour_iter_num = interface.current_iteration[neighbour]
+
+            Lambda = self.lambda_param
+            if iteration == 0:
+                n = df.FacetNormal(self.mesh)
+                pw_prev_iter = self.pressure_prev_timestep
+                k = self.relative_permeability
+                S = self.saturation
+                mu = self.viscosity
+
+                if not neighbour_iter_num == 0:
+                    # save the gl term from the neighbour first to use it in the
+                    # next iteration
+                    if self.isRichards:
+                        interface.gli_term_prev[subdomain].update(#
+                            {'wetting': interface.gli_term[neighbour]['wetting']}
+                            )
+                    else:
+                        interface.gli_term_prev[subdomain].update(#
+                            {'wetting': interface.gli_term[neighbour]['wetting']},#
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+                            )
+
+                if self.isRichards:
+                    mu_w = mu['wetting']
+                    kw = k['wetting']
+                    pwn_prev = pw_prev_iter['wetting']
+                    #calc gl0 from the flux term
+                    flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
+                    gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
+                    interface.gli_term[subdomain].update({'wetting': gl0})
+                else:
+                    mu_w = mu['wetting']
+                    mu_nw = mu['nonwetting']
+                    pwn_prev = pw_prev_iter['wetting']
+                    pnwn_prev = pw_prev_iter['nonwetting']
+                    kw = k['wetting']
+                    knw = k['nonwetting']
+                    #calc gl0 from the flux term
+                    fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
+                    fluxnw = -1/mu_nw*knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
+                    gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
+                    gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
+                    interface.gli_term[subdomain].update({'wetting': gwl0})
+                    interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+
+                interface.current_iteration[subdomain] += 1
+
+            # else: # iteration > 0
+            #     if neighbour_iter_num == iteration:
+            #         if self.isRichards:
+            #             interface.gli_term_prev[subdomain].update(#
+            #                 {'wetting': interface.gli_term[neighbour]['wetting']}
+            #                 )
+            #
+            #             pwn_prev = pw_prev_iter['wetting']
+            #             #calc gl0 from the flux term
+            #             flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
+            #             gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
+            #             interface.gli_term[subdomain].update({'wetting': gl0})
+            #         else:
+            #             mu_w = mu['wetting']
+            #             mu_nw = mu['nonwetting']
+            #             pwn_prev = pw_prev_iter['wetting']
+            #             pnwn_prev = pw_prev_iter['nonwetting']
+            #             kw = k['wetting']
+            #             knw = k['nonwetting']
+            #             #calc gl0 from the flux term
+            #             fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
+            #             fluxnw = -1/mu_nw*knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
+            #             gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
+            #             gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
+            #             interface.gli_term[subdomain].update({'wetting': gwl0})
+            #             interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+            #
+            #     else:
+            #         print('uiae')
+            #
+            #     interface.current_iteration[subdomain] += 1
 
 # END OF CLASS DomainPatch