diff --git a/.gitignore b/.gitignore
index 62fd67bc..f9d2398c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,8 +15,7 @@ __pycache__/
 *.VC.opendb
 
 # Distribution / packaging
-.Python
-build/
+.Pythonbuild/
 develop-eggs/
 dist/
 downloads/
@@ -100,3 +99,4 @@ cmake-build-*/
 *.sbn
 *.sbx
 nohup.out
+/demo/results/
diff --git a/cmf/geometry/irregular_grid.py b/cmf/geometry/irregular_grid.py
new file mode 100644
index 00000000..8f729af2
--- /dev/null
+++ b/cmf/geometry/irregular_grid.py
@@ -0,0 +1,202 @@
+import numpy as np
+import scipy.spatial as scsp
+from shapely.geometry import Polygon
+
+from .. import project as Project
+from . import create_cell, mesh_project
+
+def delaunay_cells(project: Project, x, y, z):
+    """
+    Creates a mesh using delaunay triangles between the given coordinates as cells.
+    cell.x, cell.y and cell.z represent the centroid of the cell. This results in loss of information.
+    In most cases one would use the voronoi_cells function
+
+    x, y and z must have the same length
+
+    :param project:
+    :param x: sequence of x coordinates
+    :param y: sequence of y coordinates
+    :param z: sequence of heights
+    :return: list of cells
+    """
+    points = np.array([x, y]).T
+    tri = scsp.Delaunay(points)
+
+    def get_edge_length(tri: scsp.Delaunay, tri1: int, tri2: int):
+        try:
+            n1, n2 = set(tri.simplices[tri1]) & set(tri.simplices[tri2])
+        except ValueError:
+            raise ValueError(
+                f'Triangle {tri1}:{tri.simplices[tri1]} and {tri2}:{tri.simplices[tri2]} do not share an edge')
+        c1, c2 = tri.points[[n1, n2]]
+        return np.sqrt(((c1 - c2) ** 2).sum())
+
+    def get_polygons(tri: scsp.Delaunay):
+
+        return [
+            Polygon(tri.points[s])
+            for s in tri.simplices
+        ]
+
+    def get_heights(tri: scsp.Delaunay, z):
+        return z[tri.simplices].mean(1)
+
+    polys = get_polygons(tri)
+    z_cell = get_heights(tri, z)
+
+    cells = [
+        create_cell(project, polygon, height)
+        for polygon, height in zip(polys, z_cell)
+    ]
+
+    for i, c in enumerate(cells):
+        for n in tri.neighbors[i]:
+            if n >= 0:
+                l = get_edge_length(tri, i, n)
+                c.topology.AddNeighbor(cells[n], l)
+
+    return cells
+
+
+# %%
+def __voronoi_finite_polygons_2d(vor, radius=None):
+    """
+    Reconstruct infinite voronoi regions in a 2D diagram to finite
+    regions.
+
+    Parameters
+    ----------
+    vor : Voronoi
+        Input diagram
+    radius : float, optional
+        Distance to 'points at infinity'.
+
+    Returns
+    -------
+    regions : list of tuples
+        Indices of vertices in each revised Voronoi regions.
+    vertices : list of tuples
+        Coordinates for revised Voronoi vertices. Same as coordinates
+        of input vertices, with 'points at infinity' appended to the
+        end.
+
+    """
+
+    if vor.points.shape[1] != 2:
+        raise ValueError("Requires 2D input")
+
+    new_regions = []
+    new_vertices = vor.vertices.tolist()
+
+    center = vor.points.mean(axis=0)
+    if radius is None:
+        radius = vor.points.ptp(axis=0).max()
+
+    # Construct a map containing all ridges for a given point
+    all_ridges = {}
+    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
+        all_ridges.setdefault(p1, []).append((p2, v1, v2))
+        all_ridges.setdefault(p2, []).append((p1, v1, v2))
+
+    # Reconstruct infinite regions
+    for p1, region in enumerate(vor.point_region):
+        vertices = vor.regions[region]
+
+        if all(v >= 0 for v in vertices):
+            # finite region
+            new_regions.append(vertices)
+            continue
+        elif p1 not in all_ridges:
+            continue
+        # reconstruct a non-finite region
+        ridges = all_ridges[p1]
+        new_region = [v for v in vertices if v >= 0]
+
+        for p2, v1, v2 in ridges:
+            if v2 < 0:
+                v1, v2 = v2, v1
+            if v1 >= 0:
+                # finite ridge: already in the region
+                continue
+
+            # Compute the missing endpoint of an infinite ridge
+
+            t = vor.points[p2] - vor.points[p1]  # tangent
+            t /= np.linalg.norm(t)
+            n = np.array([-t[1], t[0]])  # normal
+
+            midpoint = vor.points[[p1, p2]].mean(axis=0)
+            direction = np.sign(np.dot(midpoint - center, n)) * n
+            far_point = vor.vertices[v2] + direction * radius
+
+            new_region.append(len(new_vertices))
+            new_vertices.append(far_point.tolist())
+
+        # sort region counterclockwise
+        vs = np.asarray([new_vertices[v] for v in new_region])
+        c = vs.mean(axis=0)
+        angles = np.arctan2(vs[:, 1] - c[1], vs[:, 0] - c[0])
+        new_region = np.array(new_region)[np.argsort(angles)]
+
+        # finish
+        new_regions.append(new_region.tolist())
+
+    return new_regions, np.asarray(new_vertices)
+
+def voronoi_polygons(x, y, buffer=0.0, mask=None):
+    """
+    Creates closed voronoi polygons for x and y coordinates. With voronoi_cells, cmf cells are easily created
+    from x, y and z coordinates. Use this function directly only, if you need to process the voronoi polygons further,
+    eg. with unions of landuse or soil maps.
+
+    If this throws an unclear error, make sure you do not have any duplicate coordinates
+
+    :param x: a sequence of x values
+    :param y: a sequence of y values
+    :param buffer: a float indicating the buffer around the hull. Expressed as a fraction of the sqrt of the hull's area
+    :param mask: a polygon to use as a mask. If given, the convex hull is ignored
+    :return: A list of polygons
+    """
+    points = np.array([x, y]).T
+    vor = scsp.Voronoi(points)
+    if not mask:
+        hull = Polygon(points[scsp.ConvexHull(points).vertices])
+        mask = hull.buffer(buffer * hull.area ** 0.5)
+    regions, vertices = __voronoi_finite_polygons_2d(vor)
+    return [Polygon(vertices[r]).intersection(mask) for r in regions]
+
+
+def voronoi_cells(project: Project, x, y, z, buffer=0.0, mask=None):
+    """
+    Creates for each (x, y, z) coordinate a voronoi (nearest neighbor) cell.
+    The cells are cut of at a buffer around the convex hull of the point or using a mask.
+
+    For each coordinate a cell is created.
+
+    If this throws an unclear error, make sure you do not have any duplicate coordinates
+
+    :param project: The cmf project
+    :param x: a sequence of x values
+    :param y: a sequence of x values
+    :param z: a sequence of x values
+    :param buffer: a float indicating the buffer around the hull. Expressed as a fraction of the sqrt of the hull's area
+    :param mask: a polygon to use as a mask. If given, the convex hull is ignored
+    :return: list of cells
+    """
+    def make_cell(project, x, y, z, shape, Id=None, with_surfacewater=False):
+        cell = project.NewCell(x, y, z, shape.area, with_surfacewater=with_surfacewater)
+        cell.geometry = shape
+        if Id is not None:
+            cell.Id = Id
+        return cell
+
+    vpolys = voronoi_polygons(x, y, z, buffer, mask)
+    cells = [
+        make_cell(project, x[i], y[i], z[i], vpolys[i], Id=i, with_surfacewater=True)
+        for i in range(len(x))
+    ]
+
+    mesh_project(project)
+
+    return cells
+
diff --git a/demo/jacobian.py b/demo/jacobian.py
new file mode 100644
index 00000000..c886f582
--- /dev/null
+++ b/demo/jacobian.py
@@ -0,0 +1,38 @@
+import cmf
+from matplotlib import pyplot as plt
+
+p = cmf.project('X')
+X, = p.solutes
+
+V1 = p.NewStorage('V1', 0, 0, 0)
+V2 = p.NewStorage('V2', 1, 0, 0)
+V3 = p.NewStorage('V3', 0, 0, 1)
+V4 = p.NewStorage('V4', 1, 0, 1)
+
+NB = p.NewNeumannBoundary('NB', V1)
+NB.flux = 1
+NB.concentration[X] = 1
+
+DB = p.NewOutlet('DB', -1, 0, 0)
+
+q_nb_1 = NB.connection_to(V1)
+q_1_2 = cmf.LinearStorageConnection(V1, V2, 0.1)
+q_1_3 = cmf.LinearStorageConnection(V1, V3, 0.5)
+q_2_4 = cmf.LinearStorageConnection(V2, V4, 1)
+q_4_3 = cmf.LinearStorageConnection(V4, V3, 0.5)
+q_3_out = cmf.LinearStorageConnection(V3, DB, 1)
+
+V2[X].decay = 0.5
+V4[X].decay = 0.5
+V3[X].decay = 0.1
+
+for S in [V1, V2, V3, V4]:
+    S.volume = 0.0
+
+solver = cmf.CVodeDiag(p)
+solver.t = cmf.Time()
+solver.integrate_until(cmf.year)
+print(solver.get_info())
+t = solver.t
+plt.imshow(solver.get_jacobian(), vmin=-1, vmax=1, cmap=plt.cm.coolwarm)
+plt.colorbar()
\ No newline at end of file
diff --git a/demo/pump_trigger.py b/demo/pump_trigger.py
index 694e762e..d21b3c22 100644
--- a/demo/pump_trigger.py
+++ b/demo/pump_trigger.py
@@ -10,9 +10,9 @@ class SchmittTriggerSourceConnection:
     It is similar to the SchmittTriggerTargetConnection,
     but that is triggered by the state of the target.
 
-    If the volume of the water source is greater then the upper regulation level,
+    If th source is greater then the upper regulation level,
     the flux is switched on, if the level is smaller then the lower
-    regulation level, it is switched off. Between those levels, the
+    regulation level, it is swie volume of the watertched off. Between those levels, the
     flux stays in the same state as before.
 
     A trigger where the flow is turned on if the source has more than
diff --git a/demo/adsorption.py b/demo/subsurface/adsorption.py
similarity index 100%
rename from demo/adsorption.py
rename to demo/subsurface/adsorption.py
diff --git a/demo/adsorption_richards.py b/demo/subsurface/adsorption_richards.py
similarity index 100%
rename from demo/adsorption_richards.py
rename to demo/subsurface/adsorption_richards.py
diff --git a/demo/cmf1d_R.py b/demo/subsurface/cmf1d_R.py
similarity index 100%
rename from demo/cmf1d_R.py
rename to demo/subsurface/cmf1d_R.py
diff --git a/demo/cmf2d.py b/demo/subsurface/cmf2d.py
similarity index 100%
rename from demo/cmf2d.py
rename to demo/subsurface/cmf2d.py
diff --git a/demo/cmf3d-shapefile.py b/demo/subsurface/cmf3d-shapefile.py
similarity index 100%
rename from demo/cmf3d-shapefile.py
rename to demo/subsurface/cmf3d-shapefile.py
diff --git a/demo/darcy_plain.py b/demo/subsurface/darcy_plain.py
similarity index 100%
rename from demo/darcy_plain.py
rename to demo/subsurface/darcy_plain.py
diff --git a/demo/macropore.py b/demo/subsurface/macropore.py
similarity index 100%
rename from demo/macropore.py
rename to demo/subsurface/macropore.py
diff --git a/demo/simple_squared_3d.py b/demo/subsurface/simple_squared_3d.py
old mode 100755
new mode 100644
similarity index 100%
rename from demo/simple_squared_3d.py
rename to demo/subsurface/simple_squared_3d.py
diff --git a/demo/simpleinfiltration.py b/demo/subsurface/simpleinfiltration.py
similarity index 100%
rename from demo/simpleinfiltration.py
rename to demo/subsurface/simpleinfiltration.py
diff --git a/demo/tracer1d.py b/demo/subsurface/tracer1d.py
similarity index 100%
rename from demo/tracer1d.py
rename to demo/subsurface/tracer1d.py
diff --git a/demo/surface/2d-hillslope-surface-runoff.py b/demo/surface/2d-hillslope-surface-runoff.py
new file mode 100644
index 00000000..90b15359
--- /dev/null
+++ b/demo/surface/2d-hillslope-surface-runoff.py
@@ -0,0 +1,243 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Dec  9 10:11:47 2022
+
+@author: philipp
+"""
+
+import cmf
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.animation import FuncAnimation
+import time
+
+cmf.set_parallel_threads(2)
+l_per_s = 86.4
+
+
+def random_grid(length=6, width=3, cellsize=0.1, slope=0.1, roughness=0.01, gap_fraction=0.0):
+    """
+    Creates a random grid
+
+    Parameters
+    ----------
+    length : float, optional
+        Length of slope in m. The default is 6.
+    width : float, optional
+        width of slope in m. The default is 3.
+    cellsize : float, optional
+        Resolution in m. The default is 0.1.
+    slope : float, optional
+        Slope of the hillslope in m/m. The default is 0.1.
+    roughness : float, optional
+        Standard deviation of noise in m. The default is 0.01.
+    gap_fraction : float, optional
+        The frequency of obstacles 0..1 [-]. The default is 0 - no obstacles.
+
+    Returns
+    -------
+    np.ndarray
+        An 3d array[3,...,...] with the x, y and z 2d grids.
+
+    """
+    y, x = np.mgrid[0:length:cellsize, 0:width:cellsize]
+    z = y * -slope + np.random.normal(scale = roughness, size=x.shape)
+    if gap_fraction:
+        gaps = np.random.uniform(size=x.shape) < gap_fraction
+        z[gaps] = np.NaN
+    
+    return np.array([x, y, z])
+    
+
+class IWGBerghausen:
+    
+    def __init__(self, grid, runoff=23.0, solver_class=None, solute=None):
+        
+        self.project = cmf.project(solute or '')
+        
+        self.cells, self.shape = self.create_grid(grid) 
+        
+        self.create_upper_boundary(runoff)
+        self.create_lower_boundary()
+        solver_class = solver_class or cmf.CVodeAdams
+
+        self.solver = solver_class(self.project)
+        if hasattr(self.solver, 'initialize'):
+            self.solver.initialize()
+        
+    def scaling(self):
+        cellsize = self.cells[0,1].x - self.cells[0,0].x, self.cells[1,0].y - self.cells[0,0].y
+        return cellsize + (
+            self.cells[0, self.shape[1]-1].x - self.cells[0, 0].x + cellsize[0],
+            self.cells[self.shape[0]-1, 0].y - self.cells[0, 0].y + cellsize[1]
+        )
+        
+        
+    def create_grid(self, grid):
+        x, y, z = grid[:3]
+        cellsize_x = x[0, 1] - x[0, 0]
+        cellsize_y = y[1, 0] - y[0, 0]
+        cells = {
+            (r, c): self.project.NewCell(x[r, c], y[r, c], z[r, c], cellsize_x * cellsize_y, True)
+            for c in range(y.shape[1])
+            for r in range(y.shape[0])
+            if not np.isnan(z[r, c])
+        }
+        
+        for r, c in cells:
+            for pos, cellsize in zip([(r, c-1), (r-1,c)], [cellsize_x, cellsize_y]):
+                if pos in cells:
+                    cells[r, c].topology.AddNeighbor(cells[pos], cellsize)
+        
+        if grid.shape[0] > 3:
+            for r,c in cells:
+                cells[r, c].surfacewater.depth = grid[-1, r, c]
+                
+        for X in self.project.solutes:
+            for r, c in cells:
+                cells[r, c].surfacewater[X].conc = 0.1
+                cells[r, c].surfacewater[X].set_abs_errtol(100)
+            
+            
+        cmf.connect_cells_with_flux(self.project.cells, cmf.DiffusiveSurfaceRunoff)
+        for c in self.project.cells:
+            c.surfacewater.puddledepth = 0.0
+            c.surfacewater.nManning = 0.05
+                    
+        return cells, x.shape
+    
+    def create_upper_boundary(self, runoff):
+        cellsize_x, cellsize_y, width, length = self.scaling()
+        # Create Neumann boundary in center of first row
+        self.upper_boundary = self.project.NewNeumannBoundary('Pump', self.project.cells[0].surfacewater) 
+        self.upper_boundary.position = (width / 2, - cellsize_x / 2, 0)
+        self.upper_boundary_cells = [self.cells[r, c] for r, c in self.cells if r==0 ]
+        # Transform runoff in l/s total to m³/day/cell
+        self.upper_boundary.flux = runoff * l_per_s / len(self.upper_boundary_cells)
+        for X in self.project.solutes:
+            self.upper_boundary.concentration[X] = 1.0
+            
+        for c in self.upper_boundary_cells:
+            self.upper_boundary.connect_to(c.surfacewater)
+            
+    def create_lower_boundary(self):
+        cellsize_x, cellsize_y, width, length = self.scaling()
+       
+        self.lower_boundary = self.project.NewOutlet(
+            'Outlet', 
+            width / 2, #  x
+            length,    #  y
+            min(c.z for c in self.project)  # z
+        )
+        self.lower_boundary_cells = [self.cells[r, c] for r, c in self.cells if r == self.shape[0] - 1 ]
+        
+        for c in self.lower_boundary_cells:
+            cmf.WaterbalanceFlux(c.surfacewater, self.lower_boundary)
+            
+    def getoutflow(self):
+        return self.lower_boundary(self.solver.t) / l_per_s
+
+        
+    def getdepth(self):
+        res = np.zeros(self.shape, dtype=float) * np.NAN
+        for pos, cell in self.cells.items():
+            res[pos] = cell.surfacewater.depth
+        return res
+    
+    def getconc(self):
+        X = self.project.solutes[0]
+        res = np.zeros(self.shape, dtype=float) * np.NAN
+        for pos, cell in self.cells.items():
+            res[pos] = cell.surfacewater[X].conc
+        return res
+        
+            
+    def gethead(self):
+        res = np.zeros(self.shape, dtype=float) * np.NAN
+        for pos in self.cells:
+            cell = self.cells[pos]
+            res[pos] = cell.surfacewater.depth + cell.z
+        return res
+    
+    def getflux(self, t=None):
+        return cmf.cell_flux_directions(self.project, t or self.solver.t)
+    
+    def extent(self):
+        cs05 = self.cellsize / 2
+        return -cs05, self.width - cs05, -cs05, self.length - cs05
+    
+class Plot:
+    def __init__(self, model, callback=None, 
+                 clim=None, cmap=plt.cm.Blues,
+                 scale=None):
+        self.model = model
+        self.callback = callback or self.model.getdepth
+        p = model.project
+        self.fig, self.ax = plt.subplots(figsize=(8,12))
+        data = self.callback()
+        self.clim = clim or (np.nanmin(data), np.nanmax(data))
+        self.img = self.ax.imshow(data, cmap=cmap, aspect='equal', origin='lower',
+                        vmin=self.clim[0], vmax=self.clim[1], interpolation='none',
+                        extent=model.extent())
+        pos = cmf.cell_positions(p.cells)
+        fluxdir = model.getflux()
+        self.quiver = self.ax.quiver(pos.X, pos.Y,
+                           fluxdir.X, fluxdir.Y,
+                           scale=scale, zorder=10, pivot='middle', color='w')
+        self.title = self.ax.set_title(model.solver.t, size=16)
+        self.colorbar = self.fig.colorbar(self.img)
+        
+    def __call__(self, t=None):        
+        self.img.set_array(self.callback())
+        self.colorbar.vmin = self.clim[0]
+        self.colorbar.vmax = self.clim[1]
+        self.img.set_clim(self.clim)
+        t = t or self.model.solver.t
+        fluxdir = cmf.cell_flux_directions(self.model.project, t)
+        self.quiver.set_UVC(fluxdir.X, fluxdir.Y)
+        self.title.set_text('%s/%s %g l/s' % (t, self.model.solver.dt, self.model.getoutflow()))
+        print(self.title.get_text())
+        return self.img, self.quiver, self.title
+    
+    def animate(self, until, dt=cmf.sec):
+        
+        iterator = iter(self.model.solver.run(cmf.Time(), until, dt))
+        return FuncAnimation(self.fig, self, frames=iterator, repeat=False)
+
+
+def run_to_array(m, until, dt):
+    t0 = time.time()
+    res = np.zeros((int(until / dt), *m.shape)) * np.NAN
+    def getvalue(m):
+        if m.project.solutes:
+            return m.getconc()
+        else:
+            return m.getdepth()
+    res[0] = getvalue(m)
+    for i, t in enumerate(
+            m.solver.run(cmf.sec * 0, until, dt)
+        ):
+        print(f'{t} : out: {m.getoutflow():0.3f}')
+        res[i] = getvalue(m)
+    print('-' * 40)
+    tt = time.time() - t0
+    print(f'{tt:0.2f} s for {m.shape} cells')
+    return res
+        
+        
+if __name__ == '__main__':
+
+    until = cmf.sec * 60
+    dt = cmf.sec * 0.1
+    #grid = random_grid(6, 3, 0.1, slope=0.1, roughness=0.01, gap_fraction=0)
+    grid = np.load('berghausen_25cm_steadystate.npy')
+    m = IWGBerghausen(grid, runoff=10, solver_class=cmf.CVodeKLU, solute='X')
+    #a = run_to_array(m, until, dt)
+    print(type(m.solver).__name__)
+    print(m.solver.get_info())
+    # p = Plot(m, clim=(0,0.1), cmap=plt.cm.cividis, scale=3000)
+
+    # anim = p.animate(cmf.sec * 60, cmf.sec * 0.1)
+    # anim.save('waterdepth.apng')
+    #plt.show()
+    # np.save('waterdepth.npy', run_to_array(until, dt))
\ No newline at end of file
diff --git a/demo/abstract_reach_network.py b/demo/surface/abstract_reach_network.py
similarity index 96%
rename from demo/abstract_reach_network.py
rename to demo/surface/abstract_reach_network.py
index 6227f4d2..1f787cfd 100644
--- a/demo/abstract_reach_network.py
+++ b/demo/surface/abstract_reach_network.py
@@ -3,7 +3,7 @@
 import time
 
 def SoluteWaterSolver(states, epsilon=1e-9):
-    wsolver = cmf.CVodeKrylov(epsilon)
+    wsolver = cmf.CVodeKLU(epsilon)
     ssolver = cmf.CVodeAdams(epsilon)
     return cmf.SoluteWaterIntegrator(states.solutes, wsolver, ssolver, states)
 
@@ -81,7 +81,7 @@ def outflux(self):
     from pylab import imshow, show
 
     print(f'{"solver":<25}level size {"init sec":<10}{"run sec":<10}{"method calls":<15}')
-    for solver_type in [cmf.CVodeKLU, cmf.CVodeKrylov, cmf.CVodeAdams]:
+    for solver_type in [SoluteWaterSolver, cmf.CVodeKLU, cmf.CVodeKrylov, cmf.CVodeAdams]:
         for level in range(1, 12):
             tstart = time.time()
             model = Model(level, solver_type, 100)
diff --git a/demo/nash_cascade.py b/demo/surface/nash_cascade.py
similarity index 100%
rename from demo/nash_cascade.py
rename to demo/surface/nash_cascade.py
diff --git a/demo/parkinglot.py b/demo/surface/parkinglot.py
similarity index 100%
rename from demo/parkinglot.py
rename to demo/surface/parkinglot.py
diff --git a/demo/surface/single_channel.py b/demo/surface/single_channel.py
new file mode 100644
index 00000000..7267f859
--- /dev/null
+++ b/demo/surface/single_channel.py
@@ -0,0 +1,211 @@
+import json
+import time
+
+import pandas as pd
+
+import cmf
+import typing
+
+class River(cmf.project):
+    @property
+    def slope(self):
+        """
+        Gets the mean slope of the river
+        """
+        return (self.reaches[0].position.z - self.reaches[-1].position.z) / (self.reaches[0].position.x - self.reaches[-1].position.x)
+
+    @property
+    def length(self):
+        """
+        Get the total length of the river in m
+        """
+        return sum(r.length for r in self.reaches)
+
+    @property
+    def depth(self):
+        """
+        Get or set the depth of the river for each reach. Set to a scalar
+        will make the depth equal for each reach, or set a sequence with a value for each reach
+        """
+        return [r.depth for r in self.reaches]
+
+    @depth.setter
+    def depth(self, value):
+        try:
+            len(value)
+        except TypeError:
+            value = [value] * len(self.reaches)
+        for r, v in zip(self.reaches, value):
+            r.depth = v
+
+    def velocity(self, t: cmf.Time) -> typing.List[float]:
+        """
+        Calculate the velocity for each reach in m/s
+        """
+        area = lambda r: r.channel.get_flux_crossection(r.depth)
+        return [
+            up.flux_to(down, t) / area(up) / 86_400 if area(up)>0 else 0.0
+            for up, down in zip(self.reaches[:-1], self.reaches[1:])
+        ]
+
+    def _make_reach(self, i, slope, reachlen):
+        """
+        Creates a reach at position i
+        """
+        channel = cmf.TriangularReach(reachlen)
+        r = self.NewReach(i * reachlen, 0, - i * reachlen * slope, channel)
+        r.Name = f'R{i:03}'
+        return r
+
+    def _make_boundaries(self, upper_boundary_flow):
+        """
+        Creates the upper and lower boundary condition
+        """
+        upper_bc = self.NewNeumannBoundary('inflow', self.reaches[0])
+        upper_bc.flux = upper_boundary_flow
+        out_x = sum(r.length for r in self.reaches)
+        lower_bc = self.NewOutlet('outflow', out_x, 0, out_x * self.slope)
+        return upper_bc, lower_bc
+
+    def _connect(self, diffusive):
+        """
+        Connects the reaches and the lower boundary condition
+        """
+        for up, down in zip(self.reaches[:-1], self.reaches[1:]):
+            up.set_downstream(down)
+        down: cmf.Reach = self.reaches[-1]
+        cmf.Manning_Kinematic(down, self.lower_bc, down.channel)
+
+    def get_steadystate_depth(self, index=0):
+        """
+        Calculates the steady state depth for the reach at index according to its influx.
+        Uses scipy's root finding function brentq
+        """
+        from scipy.optimize import brentq
+        r = self.reaches[index]
+        def f(d):
+            r.depth = d
+            return r(cmf.Time())
+        return brentq(f, 0, 100)
+
+    def __init__(self, n_reaches, totlength, slope, upper_boundary_flow, solutes='', diffusive=False):
+        """
+        Creates a project for a river with constant slope, no forks and only inflow at the upper end
+        :param n_reaches: Number of numerical reaches
+        :param totlength: Total length of the river in m
+        :param slope: Constant slope of the river in m/m
+        :param upper_boundary_flow: in flow from the upper boundary in m³/day
+        :param solutes: A string denoting the names of solutes, if any
+        :param diffusive: If True, use cmf.Manning_Diffusive for flow between the reaches. Outflow to boundary is always
+        Kinematic
+        """
+        super().__init__(solutes)
+        self.reachlen = totlength / n_reaches
+        reaches = [
+            self._make_reach(i, slope, self.reachlen)
+            for i in range(n_reaches)
+        ]
+        self.upper_bc, self.lower_bc = self._make_boundaries(upper_boundary_flow)
+        self._connect(diffusive)
+
+def run_conc_step_function(p: River, solver, verbose=True, concentration=1):
+    """
+    Runs a concentration step function. River should be in steady state for water and
+    the solute storage should be zero for all storages.
+
+    First, the theoretical time for the solute to run to the outlet is calculated:
+
+    .. math ::
+        \\tau = \\frac{l}{v}
+
+    The upper boundary condition's concentration is set to `concentration` (usually 1, except for scaling experiments).
+    The model is run until :math:`2 \\tau` and saves the time, when [X]>p*C_0 as t_p, where p is [.01, .50, .99]
+    If the level is never reach t_p is 2*tau. The difference
+    between tau (t_theory), t_01 and t_99 is an indicator for numerical dispersion
+
+    :param p: The river project
+    :param solver: A solver
+    :param verbose: If True, print to stdout during the model runtime
+    :return: t_theory, t_01, t_50, t_99
+    """
+    X = p.solutes[0]
+    p.upper_bc.concentration[X] = concentration
+    def max_reach(threshold):
+        return max([i for i, r in enumerate(p.reaches) if r[X].conc > threshold] or [-1])
+    t_theory = cmf.sec * (p.length / p.velocity(solver.t)[-1])
+    t_01 = t_50 = t_99 = cmf.Time()
+    for t in solver.run(cmf.Time(), t_theory * 2, cmf.min):
+        lbcX = p.lower_bc.conc(t, X)
+        if lbcX < 0.01*concentration: t_01 = t
+        if lbcX < 0.50*concentration: t_50 = t
+        if lbcX < 0.99*concentration: t_99 = t
+        if verbose:
+            print(f'{t:%H:%M:%S} : R{max_reach(.01)}[X]>.01, R{max_reach(.99)}[X]>.99')
+    return t_theory, t_01, t_50, t_99
+
+
+class Result(dict):
+    """
+    Creates a result dictionary from the current state of the solver
+    :param solver:
+    :param kwargs:
+    :return:
+    """
+    def __init__(self, p: cmf.project, solver: cmf.CVodeBase, **kwargs):
+        t = solver.t
+        info = solver.get_info()
+        self.update(
+            dict(
+                n_reaches=len(p.reaches),
+                length=p.length,
+                slope=p.slope,
+                inflow=p.upper_bc(t),
+                solutes= ', '.join(f'{s}' for s in p.solutes),
+                t=t / cmf.sec,
+                solver=solver.to_string(),
+                open_mp_threads=cmf.get_parallel_threads()
+            ) | {
+                k: getattr(info, k)
+                for k in dir(info)
+                if not any(k.startswith(pre) for pre in ('this', 'to_', '_', 'sundials'))
+            } | kwargs
+        )
+
+    def __str__(self):
+        return '\n'.join(f'{k}: {v}' for k, v in result.items())
+
+
+def main(n_reaches, totlength, slope, upper_boundary_flow, solutes='', diffusive=False, verbose=False):
+    # Make model
+    p = River(n_reaches, totlength, slope, upper_boundary_flow, solutes=solutes, diffusive=diffusive)
+    # Bring in steady state (with empty solute storages)
+    p.depth = p.get_steadystate_depth()
+    # Create frst solver
+    solver = cmf.CVodeKLU(p)
+    # Run to ensure the steady state
+    solver(cmf.h)
+    wallclock_start = time.time()
+    cpu_start = time.process_time()
+    # Make a second solver
+    solver = cmf.CVodeKLU(p)
+    # Run the step function
+    t_theory, t_01, t_50, t_99 = run_conc_step_function(p, solver, verbose=verbose)
+    return Result(
+        p, solver,
+        step_t_theory=t_theory/cmf.sec, step_t_01=t_01/cmf.sec, step_t_50=t_50/cmf.sec, step_t_99=t_99/cmf.sec,
+        wallclock_time=time.process_time() - wallclock_start,
+        cpu_time=time.process_time() - cpu_start
+    )
+
+
+
+if __name__ == '__main__':
+    l_per_s = 86.4  #  m³/day
+    cmf.set_parallel_threads(1)
+
+    for i in range(1, 10):
+        result = main(2 ** i, 1024, 0.001, 10 * l_per_s, solutes='X')
+        print(2 ** i)
+        pd.DataFrame([result], index=[i]).to_hdf('data/single_channel.h5', 'table', append=True)
+    df = pd.read_hdf('data/single_channel.h5', 'table')
+    print(df)
\ No newline at end of file
diff --git a/demo/surface/voronoi_demo.py b/demo/surface/voronoi_demo.py
new file mode 100644
index 00000000..4db1ea6c
--- /dev/null
+++ b/demo/surface/voronoi_demo.py
@@ -0,0 +1,49 @@
+import cmf
+import cmf.geometry
+from cmf.geometry import irregular_grid as cmf_ig
+
+import numpy as np
+from shapely.geometry import Polygon
+
+def create_points(n_cells, width, height, slope, roughness):
+    x = np.random.uniform(0, width, n_cells)
+    y = np.random.uniform(0, height, n_cells)
+    z = y * slope + np.random.normal(scale=roughness, size=n_cells)
+    return x, y, z
+
+
+def run(project, until, dt, solver=cmf.CVodeKLU):
+    for c in project:
+        c.surfacewater.depth = 0.0
+    if type(solver) is type:
+        solver = solver(project)
+    for t in solver.run(solver.t, until, dt):
+        print(t)
+    print(solver.get_info())
+    return solver
+
+
+def add_boundaries(project: cmf.project, flux=10):
+    in_cell = max(project, key=lambda c: c.z)
+    out_cell = min(project, key=lambda c: c.z)
+
+    in_boundary = project.NewNeumannBoundary('inflow', in_cell.surfacewater)
+    out_boundary = project.NewOutlet('outflow', out_cell.x, out_cell.y - 0.5, out_cell.z - 0.01)
+    cmf.DiffusiveSurfaceRunoff(out_cell.surfacewater, out_boundary, 1)
+
+    return in_boundary, out_boundary
+
+
+if __name__ == '__main__':
+
+    x, y, z = create_points(250, 10, 10, 0.1, 0.01)
+    p = cmf.project()
+    cells = cmf_ig.voronoi_cells(p, x, y, z, 0.05)
+    cmf.connect_cells_with_flux(p, cmf.DiffusiveSurfaceRunoff)
+    b_in, b_out = add_boundaries(p, 10.0)
+
+    solver = run(p, cmf.min, cmf.sec, cmf.CVodeKLU)
+
+
+
+
diff --git a/demo/tracer_filter.py b/demo/tracer_filter.py
index d67d041e..36a915d0 100644
--- a/demo/tracer_filter.py
+++ b/demo/tracer_filter.py
@@ -1,11 +1,9 @@
-# -*- coding: utf-8 -*-
 """
-Created on Tue May  2 12:33:44 2017
-
-@author: gh1961
+Shows the effect of a solute filtering flux connection
 """
-# %%
 import cmf
+import datetime
+
 
 p = cmf.project('no_filter filter')
 NF, F = p.solutes
@@ -19,20 +17,20 @@
 q.set_tracer_filter(F, 0.5)
 # %%
 W1.volume = 1.0
-W1.conc(NF, 1.0)
-W1.conc(F, 1.0)
+W1[NF].state = 1.0
+W1[F].state = 1.0
+
 W2.volume = 0.0
-W2.Solute(NF).state = 0.0
-W2.Solute(F).state = 0.0
+W2[NF].state = 0.0
+W2[F].state = 0.0
 # Create an integrator for the ODE represented by project p, with an error tolerance of 1e-9
 solver = cmf.CVodeAdams(p, 1e-9)
 # Import Python's datetime module
-import datetime
 
 # Set the intitial time of the solver
 solver.t = datetime.datetime(2012, 1, 1)
 # %%
-result = [[W1.Solute(NF).state, W2.Solute(NF).state, W1.Solute(F).state, W2.Solute(F).state] for t in
+result = [[W1[NF].state, W2[NF].state, W1[F].state, W2[F].state] for t in
           solver.run(datetime.datetime(2012, 1, 1), datetime.datetime(2012, 1, 7), datetime.timedelta(hours=1))]
 import pylab as plt
 
diff --git a/demo/CmfTutChannel.py b/demo/tutorial/CmfTutChannel.py
similarity index 100%
rename from demo/CmfTutChannel.py
rename to demo/tutorial/CmfTutChannel.py
diff --git a/demo/CmfTutET.py b/demo/tutorial/CmfTutET.py
similarity index 100%
rename from demo/CmfTutET.py
rename to demo/tutorial/CmfTutET.py
diff --git a/demo/cmfTutBoundary.py b/demo/tutorial/cmfTutBoundary.py
similarity index 100%
rename from demo/cmfTutBoundary.py
rename to demo/tutorial/cmfTutBoundary.py
diff --git a/demo/lumped_2storages.py b/demo/tutorial/lumped_2storages.py
similarity index 100%
rename from demo/lumped_2storages.py
rename to demo/tutorial/lumped_2storages.py
diff --git a/demo/simple_connections.py b/demo/tutorial/simple_connections.py
similarity index 100%
rename from demo/simple_connections.py
rename to demo/tutorial/simple_connections.py