diff --git a/LDDsimulation.py b/LDDsimulation.py
index b2e2d112a8854a0a62a00786d7f9963bdc92d471..a6848c82ef56b76ad6ec04f6106132f5e31e3b11 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -86,6 +86,8 @@ class LDDsimulation(object):
         self._parameters_set = False
         # dictionary of objects of class DomainPatch initialised by self._init_subdomains()
         self.subdomain = dict()
+        # dictionary to store the dirichletBC in.
+        self.outerBC = dict()
 
     def set_parameters(self,#
                        output_dir: str,#
@@ -103,7 +105,8 @@ class LDDsimulation(object):
                        saturation: tp.Dict[int, tp.Callable[...,None]],#
                        timestep: tp.Dict[int, float],#
                        sources: tp.Dict[int, tp.Dict[str, str]],#
-                       initial_conditions: tp.Dict[int, tp.Dict[str, str]]
+                       initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
+                       outer_dirichletBC: tp.Dict[int, tp.Dict[str, str]],#
                        )-> None:
 
         """ set parameters of an instance of class LDDsimulation"""
@@ -122,6 +125,7 @@ class LDDsimulation(object):
         self.timestep = timestep
         self.sources = sources
         self.initial_conditions = initial_conditions
+        self.outer_dirichletBC = outer_dirichletBC
         self._parameters_set = True
 
     def initialise(self) -> None:
@@ -136,6 +140,8 @@ class LDDsimulation(object):
         self._init_subdomains()
         self._init_initial_values()
         self._eval_sources()
+        # init Dirichlet Boundary values for the first time
+        self.set_DirichletBC(time = 0, first_init = True)
 
     ## Private methods
     def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
@@ -328,26 +334,32 @@ class LDDsimulation(object):
                 subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
                                                     'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
 
-    def _init_DirichletBC(self, interpolation_degree: int = 2):
-        """ set initial values
+    def set_DirichletBC(self, interpolation_degree: int = 2, #
+                        time: float = 0,#
+                        # exact_solution: bool = False,#
+                        first_init: bool = False):
+        """ set dirichlet boundary values at time time.
         """
         for num, subdomain in self.subdomain.items():
             mesh = subdomain.mesh
             V = subdomain.function_space
-            # p0 contains both pw_0 and pnw_0 for subdomain[num]
-            p0 = self.initial_conditions[num]
+            # if first_init==True we have not initialised the dictionary
+            if first_init:
+                self.outerBC.update({num: dict()})
+
+            pDC = self.outer_dirichletBC[num]
+            boundary_marker = subdomain.boundary_marker
             if subdomain.isRichards:
                 # note that the default interpolation degree is 2
-                pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
-                subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
-                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting'])}
+                pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
+                self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 0)})
             else:
-                pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
-                pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
-                subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting']),#
-                                                'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
-                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
-                                                    'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
+                pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
+                pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
+                self.outerBC[num].update(#
+                    {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 0)},#
+                    {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 0)},#
+                )
 
     def _eval_sources(self, interpolation_degree: int = 2, time: float = 0):
         """ evaluate time dependent source terms or initialise them if time == 0