diff --git a/.gitignore b/.gitignore
index 2a29ad13..9eefc3c2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -77,6 +77,9 @@ ENV/
 env.bak/
 venv.bak/
 
+# development notes
+dev-notes/
+
 # Spyder project settings
 .spyderproject
 .spyproject
diff --git a/aeolis/model.py b/aeolis/model.py
index 22a68873..07f984d9 100644
--- a/aeolis/model.py
+++ b/aeolis/model.py
@@ -29,7 +29,6 @@
 
 import os
 import imp
-import sys
 import time
 import glob
 import logging
@@ -59,11 +58,18 @@
 import aeolis.fences
 import aeolis.gridparams
 
+# type hints
+from typing import Any, Union, Tuple
+from numpy import ndarray
+
 from aeolis.utils import *
 
 class StreamFormatter(logging.Formatter):
+    """A formater for log messages"""
+
+    def format(self, record) -> str:
+        """Formats log messages for consult standard output"""
 
-    def format(self, record):
         if record.levelname == 'INFO':
             return record.getMessage()
         else:
@@ -73,10 +79,10 @@ def format(self, record):
 # initialize logger
 logger = logging.getLogger(__name__)
 
-
 __version__ = ''
 __root__ = os.path.dirname(__file__)
 
+# Collect model version from VERSION file
 try:
     __version__ = open(os.path.join(__root__, 'VERSION')).read().strip()
 except:
@@ -84,37 +90,41 @@ def format(self, record):
 
 
 class ModelState(dict):
-    '''Dictionary-like object to store model state
-
+    """Dictionary-like object to store model state.
     Model state variables are mutable by default, but can be set
-    immutable. In the latter case any actions that set the immutable
+    to immutable. In the latter case any actions that set the immutable
     model state variable are ignored.
+    """
 
-    '''
-
-
-    def __init__(self, *args, **kwargs):
+    def __init__(self, *args, **kwargs) -> None:
+        """Initialize dictionary for model states"""
         self.ismutable = set()
         super(ModelState, self).__init__(*args, **kwargs)
 
 
-    def __setitem__(self, k, v):
+    def __setitem__(self, k:str, v:Any) -> None:
+        """sets new and existing mutable model variables and thier values 
+        to the model state
+        """
+        
         if k not in self.keys() or k in self.ismutable:
             super(ModelState, self).__setitem__(k, v)
             self.set_mutable(k)
 
 
-    def set_mutable(self, k):
+    def set_mutable(self, k:str) -> None:
+        """Adds mutable variable to the model state"""
         self.ismutable.add(k)
 
 
-    def set_immutable(self, k):
+    def set_immutable(self, k:str) -> None:
+        """Removes existing model varialble from the set of mutable variables"""
         if k in self.ismutable:
             self.ismutable.remove(k)
 
 
 class AeoLiS(IBmi):
-    '''AeoLiS model class
+    """AeoLiS model class
 
     AeoLiS is a process-based model for simulating supply-limited
     aeolian sediment transport. This model class is compatible with
@@ -137,47 +147,43 @@ class AeoLiS(IBmi):
     >>> for i in range(10):
     >>>     model.update(60.) # step 60 seconds forward
     >>> model.finalize()
+    """
 
-    '''
-
-    def __init__(self, configfile):
+    def __init__(self, configfile: str) -> None:
         '''Initialize class
 
         Parameters
         ----------
-        configfile : str
-            Model configuration file. See :func:`~inout.read_configfile()`.
+        configfile:
+            Path to model configuration file. See :func:`~inout.read_configfile()`.
 
         '''
 
-        self.t = 0.
-        self.dt = 0.
-        self.configfile = ''
+        self.t: float = 0
+        self.dt: float = 0
+        self.configfile: str = ''
 
-        self.l = {} # previous spatial grids
-        self.s = ModelState() # spatial grids
-        self.p = {} # parameters
-        self.c = {} # counters
+        self.l: dict = {} # previous spatial grids
+        self.s: ModelState = ModelState() # spatial grids
+        self.p: dict = {} # parameters
+        self.c: dict = {} # counters
 
         self.configfile = configfile
 
-
     def __enter__(self):
+        """Initialize AeoliS class"""
         self.initialize()
         return self
 
-
-    def __exit__(self, *args):
+    def __exit__(self, *args) -> None:
+        """Calls to terminate AeoLiS class"""
         self.finalize()
 
-
-    def initialize(self):
-        '''Initialize model
-
-        Read model configuration file and initialize parameters and
-        spatial grids dictionary and load bathymetry and bed
+    def initialize(self)-> None:
+        '''
+        Read model configuration file, initialize parameters and
+        spatial grids dictionary, and load bathymetry and bed
         composition.
-
         '''
 
         # read configuration file
@@ -235,7 +241,7 @@ def initialize(self):
             aeolis.inout.visualize_timeseries(self.p, self.t)
 
 
-    def update(self, dt=-1):
+    def update(self, dt:float=-1) -> None:
         '''Time stepping function
 
         Takes a single step in time. Interpolates wind and
@@ -253,7 +259,7 @@ def update(self, dt=-1):
 
         Parameters
         ----------
-        dt : float, optional
+        dt : optional
             Time step in seconds. The time step specified in the model
             configuration file is used in case dt is smaller than
             zero. For explicit numerical schemes the time step is
@@ -269,7 +275,7 @@ def update(self, dt=-1):
         self.l['dzbavg'] = self.s['dzbavg'].copy()
         
         # interpolate wind time series
-        self.s = aeolis.wind.interpolate(self.s, self.p, self.t)     
+        self.s = aeolis.wind.interpolate(self.s, self.p, self.t)  
     
         # Rotate gridparams, such that the grids is alligned horizontally
         self.s = self.grid_rotate(self.p['alpha'])
@@ -330,7 +336,8 @@ def update(self, dt=-1):
         # compute dune erosion
         if self.p['process_dune_erosion']:
             self.s = aeolis.erosion.run_ph12(self.s, self.p, self.t)
-            self.s = aeolis.avalanching.angele_of_repose(self.s, self.p) #Since the aeolian module is only run for winds above threshold, also run avalanching routine here
+            self.s = aeolis.avalanching.angele_of_repose(self.s, self.p) #Since the aeolian module is 
+            # only run for winds above threshold, also run avalanching routine here
             self.s = aeolis.avalanching.avalanche(self.s, self.p)
 
         # grow vegetation
@@ -352,47 +359,38 @@ def update(self, dt=-1):
 
     def finalize(self):
         '''Finalize model'''
-
         pass
 
 
-    def get_current_time(self):
+    def get_current_time(self) -> float:
         '''
         Returns
         -------
         float
             Current simulation time
-
         '''
 
         return self.t
 
-
-    def get_end_time(self):
+    def get_end_time(self) -> float:
         '''
         Returns
         -------
-        float
             Final simulation time
-
         '''
 
         return self.p['tstop']
 
-
-    def get_start_time(self):
+    def get_start_time(self) -> float:
         '''
         Returns
         -------
-        float
             Initial simulation time
-
         '''
 
         return self.p['tstart']
 
-
-    def get_var(self, var):
+    def get_var(self, var:str) -> Union[ndarray, int, float, str, list]:
         '''Returns spatial grid or model configuration parameter
 
         If the given variable name matches with a spatial grid, the
@@ -403,7 +401,7 @@ def get_var(self, var):
 
         Parameters
         ----------
-        var : str
+        var:
             Name of spatial grid or model configuration parameter
 
         Returns
@@ -422,7 +420,6 @@ def get_var(self, var):
         See Also
         --------
         model.AeoLiS.set_var
-
         '''
 
         if var in self.s:
@@ -436,31 +433,28 @@ def get_var(self, var):
             return None
 
 
-    def get_var_count(self):
+    def get_var_count(self) -> int:
         '''
         Returns
         -------
-        int
             Number of spatial grids
-
         '''
 
         return len(self.s)
 
 
-    def get_var_name(self, i):
+    def get_var_name(self, i: int) -> Union[str, int]:
         '''Returns name of spatial grid by index (in alphabetical order)
 
         Parameters
         ----------
-        i : int
+        i : 
             Index of spatial grid
 
         Returns
         -------
         str or -1
             Name of spatial grid or -1 in case index exceeds the number of grids
-
         '''
 
         if len(self.s) > i:
@@ -468,20 +462,17 @@ def get_var_name(self, i):
         else:
             return -1
 
-
-    def get_var_rank(self, var):
+    def get_var_rank(self, var: str) -> int:
         '''Returns rank of spatial grid
 
         Parameters
         ----------
-        var : str
+        var :
             Name of spatial grid
 
         Returns
         -------
-        int
             Rank of spatial grid or -1 if not found
-
         '''
 
         if var in self.s:
@@ -490,19 +481,17 @@ def get_var_rank(self, var):
             return -1
 
 
-    def get_var_shape(self, var):
+    def get_var_shape(self, var:str) -> Union[Tuple, str]:
         '''Returns shape of spatial grid
 
         Parameters
         ----------
-        var : str
+        var : 
             Name of spatial grid
 
         Returns
         -------
-        tuple or int
             Dimensions of spatial grid or -1 if not found
-
         '''
 
         if var in self.s:
@@ -511,19 +500,17 @@ def get_var_shape(self, var):
             return -1
 
 
-    def get_var_type(self, var):
-        '''Returns variable type of spatial grid
+    def get_var_type(self, var:str) -> Union[str, int]:
+        '''Returns data type of variable in spatial grid
 
         Parameters
         ----------
-        var : str
+        var :
             Name of spatial grid
 
         Returns
         -------
-        str or int
             Variable type of spatial grid or -1 if not found
-
         '''
 
         if var in self.s:
@@ -531,16 +518,15 @@ def get_var_type(self, var):
         else:
             return -1
 
-
-    def inq_compound(self):
+    def inq_compound(self) -> None:
         logger.log_and_raise('Method not yet implemented [inq_compound]', exc=NotImplementedError)
 
 
-    def inq_compound_field(self):
+    def inq_compound_field(self) -> None:
         logger.log_and_raise('Method not yet implemented [inq_compound_field]', exc=NotImplementedError)
 
 
-    def set_var(self, var, val):
+    def set_var(self, var:str, val:Union[ndarray, int, float, str, list]) -> None:
         '''Sets spatial grid or model configuration parameter
 
         If the given variable name matches with a spatial grid, the
@@ -550,9 +536,9 @@ def set_var(self, var, val):
 
         Parameters
         ----------
-        var : str
+        var : 
             Name of spatial grid or model configuration parameter
-        val : np.ndarray or int, float, str or list
+        val : 
             Spatial grid or model configuration parameter
 
         Examples
@@ -566,7 +552,6 @@ def set_var(self, var, val):
         See Also
         --------
         model.AeoLiS.get_var
-
         '''
 
         if var in self.s:
@@ -575,27 +560,26 @@ def set_var(self, var, val):
             self.p[var] = val
 
 
-    def set_var_index(self, i, val):
+    def set_var_index(self, i:int, val:ndarray) -> None:
         '''Set spatial grid by index (in alphabetical order)
 
         Parameters
         ----------
-        i : int
+        i :
             Index of spatial grid
-        val : np.ndarray
+        val : 
             Spatial grid
-
         '''
 
         var = self.get_var_name(i)
         self.set_var(var, val)
 
 
-    def set_var_slice(self):
+    def set_var_slice(self) -> None:
         logger.log_and_raise('Method not yet implemented [set_var_slice]', exc=NotImplementedError)
 
 
-    def set_timestep(self, dt=-1.):
+    def set_timestep(self, dt:float=-1.0) -> bool:
         '''Determine optimal time step
 
         If no time step is given the optimal time step is
@@ -619,12 +603,11 @@ def set_timestep(self, dt=-1.):
 
         Parameters
         ----------
-        df : float, optional
+        df :
             Preferred time step
 
         Returns
         -------
-        bool
             False if determination of time step was unsuccessful, True otherwise
 
         '''
@@ -655,10 +638,9 @@ def set_timestep(self, dt=-1.):
             
         self.p['dt_opt'] = self.dt
 
-
         return True
 
-    def grid_rotate(self, angle):
+    def grid_rotate(self, angle:float) -> ModelState:
         
         s = self.s
         p = self.p
@@ -680,7 +662,7 @@ def grid_rotate(self, angle):
         
         return s
     
-    def euler_forward(self):
+    def euler_forward(self) -> Any:
         '''Convenience function for explicit solver based on Euler forward scheme
 
         See Also
@@ -701,7 +683,7 @@ def euler_forward(self):
         return solve
 
 
-    def euler_backward(self):
+    def euler_backward(self) -> Any:
         '''Convenience function for implicit solver based on Euler backward scheme
 
         See Also
@@ -721,13 +703,12 @@ def euler_backward(self):
             
         return solve
 
-    def crank_nicolson(self):
+    def crank_nicolson(self) -> Any:
         '''Convenience function for semi-implicit solver based on Crank-Nicolson scheme
 
         See Also
         --------
         model.AeoLiS.solve
-
         '''
 
         if self.p['solver'].lower() == 'trunk':
@@ -742,9 +723,8 @@ def crank_nicolson(self):
         return solve
 
 
-    def solve_steadystate(self):
+    def solve_steadystate(self) -> dict:
         '''Implements the steady state solution
-
         '''
         # upwind scheme:
         beta = 1. 
@@ -993,8 +973,7 @@ def solve_steadystate(self):
                 # solve system with current weights
                 Ct_i = scipy.sparse.linalg.spsolve(A, y_i.flatten())
                 Ct_i = prevent_tiny_negatives(Ct_i, p['max_error'])
-                
-                    
+                   
                 # check for negative values
                 if Ct_i.min() < 0.:
                     ix = Ct_i < 0.
@@ -1087,7 +1066,7 @@ def solve_steadystate(self):
                     q=q)
         
         
-    def solve(self, alpha=.5, beta=1.):
+    def solve(self, alpha:float=.5, beta:float=1.) -> dict:
         '''Implements the explicit Euler forward, implicit Euler backward and semi-implicit Crank-Nicolson numerical schemes
 
         Determines weights of sediment fractions, sediment pickup and
@@ -1097,14 +1076,13 @@ def solve(self, alpha=.5, beta=1.):
 
         Parameters
         ----------
-        alpha : float, optional
+        alpha :
             Implicitness coefficient (0.0 for Euler forward, 1.0 for Euler backward or 0.5 for Crank-Nicolson, default=0.5)
-        beta : float, optional
+        beta : 
             Centralization coefficient (1.0 for upwind or 0.5 for centralized, default=1.0)
 
         Returns
         -------
-        dict
             Partial spatial grid dictionary
 
         Examples
@@ -1382,13 +1360,10 @@ def solve(self, alpha=.5, beta=1.):
                 if p['boundary_onshore'] == 'constant':
                     y_i[:,-1] = p['constant_onshore_flux'] / s['u'][:,-1,i]
 
-
-
                 # solve system with current weights
                 Ct_i = scipy.sparse.linalg.spsolve(A, y_i.flatten())
                 Ct_i = prevent_tiny_negatives(Ct_i, p['max_error'])
                 
-
                 # check for negative values
                 if Ct_i.min() < 0.:
                     ix = Ct_i < 0.
@@ -1418,7 +1393,6 @@ def solve(self, alpha=.5, beta=1.):
                 ix = (deficit_i > p['max_error']) \
                      & (w_i * Cu_i > 0.)
 
-
                 # quit the iteration if there is no deficit, otherwise
                 # back-compute the maximum weight allowed to get zero
                 # deficit for the current fraction and progress to
@@ -1470,12 +1444,9 @@ def solve(self, alpha=.5, beta=1.):
                                       minweight=np.sum(w, axis=-1).min(),
                                       **logprops))
            
-
-
         qs = Ct * s['us'] 
         qn = Ct * s['un'] 
 
-
         return dict(Ct=Ct,
                     qs=qs,
                     qn=qn,
@@ -1486,7 +1457,7 @@ def solve(self, alpha=.5, beta=1.):
                     w_bed=w_bed)
         
         
-    def solve_steadystatepieter(self):
+    def solve_steadystatepieter(self) -> dict:
         
         beta = 1. 
         
@@ -1710,8 +1681,7 @@ def solve_steadystatepieter(self):
                 yCt_i[:,:-1]  -= s['dn'][:,:-1] * ufs[:,1:-1] * Ctxfs_i[:,1:-1] #upper x-face
                 yCt_i[1:,:]   += s['ds'][1:,:]  * ufn[1:-1,:] * Ctxfn_i[1:-1,:] #lower y-face
                 yCt_i[:-1,:]  -= s['ds'][:-1,:] * ufn[1:-1,:] * Ctxfn_i[1:-1,:] #upper y-face
-             
-                    
+                
                 # boundary conditions
                 # offshore boundary (i=0)
 
@@ -1791,7 +1761,6 @@ def solve_steadystatepieter(self):
                     Ct_i[~ix] *= 1. + Ct_i[ix].sum() / Ct_i[~ix].sum()
                     Ct_i[ix] = 0.
 
-
                 # determine pickup and deficit for current fraction
                 Cu_i = s['Cu'][:,:,i].flatten()
                 mass_i = s['mass'][:,:,0,i].flatten()
@@ -1845,7 +1814,6 @@ def solve_steadystatepieter(self):
                                        **logprops))
         # end loop over frations
 
-
         # check if there are any cells where the sum of all weights is
         # smaller than unity. these cells are supply-limited for all
         # fractions. Log these events.
@@ -1869,7 +1837,7 @@ def solve_steadystatepieter(self):
                     w_bed=w_bed)
     
     
-    def solve_pieter(self, alpha=.5, beta=1.):
+    def solve_pieter(self, alpha:float=.5, beta:float=1.) -> dict:
         '''Implements the explicit Euler forward, implicit Euler backward and semi-implicit Crank-Nicolson numerical schemes
 
         Determines weights of sediment fractions, sediment pickup and
@@ -1879,14 +1847,13 @@ def solve_pieter(self, alpha=.5, beta=1.):
 
         Parameters
         ----------
-        alpha : float, optional
+        alpha : 
             Implicitness coefficient (0.0 for Euler forward, 1.0 for Euler backward or 0.5 for Crank-Nicolson, default=0.5)
         beta : float, optional
             Centralization coefficient (1.0 for upwind or 0.5 for centralized, default=1.0)
 
         Returns
         -------
-        dict
             Partial spatial grid dictionary
 
         Examples
@@ -1902,10 +1869,8 @@ def solve_pieter(self, alpha=.5, beta=1.):
         model.AeoLiS.crank_nicolson
         transport.compute_weights
         transport.renormalize_weights
-
         '''
 
-        
         l = self.l
         s = self.s
         p = self.p
@@ -2080,7 +2045,6 @@ def solve_pieter(self, alpha=.5, beta=1.):
             # solve transport for each fraction separately using latest
             # available weights
         
-
             # renormalize weights for all fractions equal or larger
             # than the current one such that the sum of all weights is
             # unity
@@ -2132,8 +2096,7 @@ def solve_pieter(self, alpha=.5, beta=1.):
                 yCt_i[:,:-1]  -= s['dn'][:,:-1] * ufs[:,1:-1] * Ctxfs_i[:,1:-1] #upper x-face
                 yCt_i[1:,:]   += s['ds'][1:,:]  * ufn[1:-1,:] * Ctxfn_i[1:-1,:] #lower y-face
                 yCt_i[:-1,:]  -= s['ds'][:-1,:] * ufn[1:-1,:] * Ctxfn_i[1:-1,:] #upper y-face
-                
-                    
+                  
                 # boundary conditions
                 # offshore boundary (i=0)
 
@@ -2215,7 +2178,6 @@ def solve_pieter(self, alpha=.5, beta=1.):
                     Ct_i[~ix] *= 1. + Ct_i[ix].sum() / Ct_i[~ix].sum()
                     Ct_i[ix] = 0.
 
-
                 # determine pickup and deficit for current fraction
                 Cu_i = s['Cu'][:,:,i].flatten()
                 mass_i = s['mass'][:,:,0,i].flatten()
@@ -2269,7 +2231,6 @@ def solve_pieter(self, alpha=.5, beta=1.):
                                        **logprops))
         # end loop over frations
 
-
         # check if there are any cells where the sum of all weights is
         # smaller than unity. these cells are supply-limited for all
         # fractions. Log these events.
@@ -2298,8 +2259,7 @@ def solve_pieter(self, alpha=.5, beta=1.):
                     w_bed=w_bed,
                     q=q)        
 
-
-    def get_count(self, name):
+    def get_count(self, name:str) -> int:
         '''Get counter value
 
         Parameters
@@ -2315,16 +2275,15 @@ def get_count(self, name):
             return 0
 
 
-    def _count(self, name, n=1):
+    def _count(self, name:str, n:int=1) -> None:
         '''Increase counter
 
         Parameters
         ----------
-        name : str
+        name : 
             Name of counter
-        n : int, optional
+        n : optional
             Increment of counter (default: 1)
-
         '''
 
         if name not in self.c:
@@ -2332,7 +2291,7 @@ def _count(self, name, n=1):
         self.c[name] += n
 
 
-    def _dims2shape(self, dims):
+    def _dims2shape(self, dims) -> Tuple:
         '''Converts named dimensions to numbered shape
 
         Supports only dimension names that can be found in the model
@@ -2347,7 +2306,6 @@ def _dims2shape(self, dims):
 
         Returns
         -------
-        tuple
             Shape of spatial grid
 
         '''
@@ -2361,12 +2319,12 @@ def _dims2shape(self, dims):
 
 
     @staticmethod
-    def dimensions(var=None):
+    def dimensions(var:str=None) -> Union[Tuple, dict]:
         '''Static method that returns named dimensions of all spatial grids
 
         Parameters
         ----------
-        var : str, optional
+        var : optional
             Name of spatial grid
 
         Returns
@@ -2376,7 +2334,6 @@ def dimensions(var=None):
             dictionary with all named dimensions of all spatial
             grids. Returns nothing if requested spatial grid is not
             defined.
-
         '''
 
         dims = {s:d
@@ -2419,11 +2376,9 @@ class to start an AeoLiS model run.
     See Also
     --------
     console.aeolis
-
     '''
 
-
-    def __init__(self, configfile='aeolis.txt'):
+    def __init__(self, configfile:str='aeolis.txt') -> None:
         '''Initialize class
 
         Reads model configuration file without parsing all referenced
@@ -2471,7 +2426,7 @@ def __init__(self, configfile='aeolis.txt'):
             logger.log_and_raise('Configuration file not found [%s]' % self.configfile, exc=IOError)
 
 
-    def run(self, callback=None, restartfile=None):
+    def run(self, callback=None, restartfile:str=None) -> None:
         '''Start model time loop
 
         Changes current working directory to the model directory,
@@ -2585,7 +2540,7 @@ def run(self, callback=None, restartfile=None):
         logging.shutdown()
 
 
-    def set_configfile(self, configfile):
+    def set_configfile(self, configfile:str) -> None:
         '''Set model configuration file name'''
 
         self.changed = False
@@ -2595,7 +2550,7 @@ def set_configfile(self, configfile):
             self.configfile = configfile
 
 
-    def set_params(self, **kwargs):
+    def set_params(self, **kwargs) -> None:
         '''Set model configuration parameters'''
 
         if len(kwargs) > 0:
@@ -2603,7 +2558,7 @@ def set_params(self, **kwargs):
             self.p.update(kwargs)
 
 
-    def get_statistic(self, var, stat='avg'):
+    def get_statistic(self, var:str, stat:str='avg') -> Union[ndarray, None]:
         '''Return statistic of spatial grid
 
         Parameters
@@ -2615,9 +2570,7 @@ def get_statistic(self, var, stat='avg'):
 
         Returns
         -------
-        numpy.ndarray
             Statistic of spatial grid
-
         '''
 
         if stat in ['min', 'max', 'sum']:
@@ -2637,7 +2590,7 @@ def get_statistic(self, var, stat='avg'):
             return None
 
 
-    def get_var(self, var, clear=True):
+    def get_var(self, var:str, clear:bool=True) -> Union[ndarray, int, float, str, list]:
         '''Returns spatial grid, statistic or model configuration parameter
 
         Overloads the :func:`~model.AeoLiS.get_var()` function and
@@ -2672,7 +2625,6 @@ def get_var(self, var, clear=True):
         See Also
         --------
         model.AeoLiS.get_var
-
         '''
 
         self.clear = clear
@@ -2693,7 +2645,7 @@ def get_var(self, var, clear=True):
         return super(AeoLiSRunner, self).get_var(var)
 
 
-    def initialize(self):
+    def initialize(self) -> None:
         '''Initialize model
 
         Overloads the :func:`~model.AeoLiS.initialize()` function, but
@@ -2704,8 +2656,7 @@ def initialize(self):
         super(AeoLiSRunner, self).initialize()
         self.output_init()
 
-
-    def update(self, dt=-1):
+    def update(self, dt:float=-1) -> None:
         '''Time stepping function
 
         Overloads the :func:`~model.AeoLiS.update()` function,
@@ -2727,7 +2678,7 @@ def update(self, dt=-1):
         self.output_update()
 
 
-    def write_params(self):
+    def write_params(self) -> None:
         '''Write updated model configuration to configuration file
 
         Creates a backup in case the model configration file already
@@ -2745,7 +2696,7 @@ def write_params(self):
             self.changed = False
 
 
-    def output_init(self):
+    def output_init(self) -> None:
 
         '''Initialize netCDF4 output file and output statistics dictionary'''
 
@@ -2781,7 +2732,7 @@ def output_init(self):
         self.output_clear()
 
 
-    def output_clear(self):
+    def output_clear(self) -> None:
         '''Clears output statistics dictionary
 
         Creates a matrix for minimum, maximum, variance and summed
@@ -2800,7 +2751,7 @@ def output_clear(self):
         self.n = 0
 
 
-    def output_update(self):
+    def output_update(self) -> None:
         '''Updates output statistics dictionary
 
         Updates matrices with minimum, maximum, variance and summed
@@ -2826,7 +2777,7 @@ def output_update(self):
         self.n += 1
 
 
-    def output_write(self):
+    def output_write(self) -> None:
         '''Appends output to netCDF4 output file
 
         If the time since the last output is equal or larger than the
@@ -2834,7 +2785,6 @@ def output_write(self):
         output file. Computes the average and variance values based on
         available output statistics and clear output statistics
         dictionary.
-
         '''
 
         if self.t - self.tout >= self.p['output_times'] or self.t == 0.:
@@ -2858,7 +2808,7 @@ def output_write(self):
             self.trestart = self.t
 
 
-    def load_hotstartfiles(self):
+    def load_hotstartfiles(self) -> None:
         '''Load model state from hotstart files
 
         Hotstart files are plain text representations of model state
@@ -2886,18 +2836,22 @@ def load_hotstartfiles(self):
                 logger.warning('Unrecognized hotstart file [%s]' % fname)
 
 
-    def load_restartfile(self, restartfile):
+    def load_restartfile(self, restartfile:str) -> bool:
         '''Load model state from restart file
 
         Parameters
         ----------
         restartfile : str
             Path to previously written restartfile.
-
+        
+        Returns
+        -------
+            True if model state from restartfile is loaded successfully
         '''
 
         if restartfile:
             if os.path.exists(restartfile):
+                # Load model state
                 with open(restartfile, 'r') as fp:
                     state = pickle.load(fp)
 
@@ -2912,12 +2866,11 @@ def load_restartfile(self, restartfile):
                     return True
             else:
                 logger.log_and_raise('Restart file not found [%s]' % restartfile, exc=IOError)
-
-        return False
+                return False
 
 
-    def dump_restartfile(self):
-        '''Dump model state to restart file'''
+    def dump_restartfile(self) -> None:
+        '''Dumps model state to restart file'''
 
         restartfile = '%s.r%010d' % (os.path.splitext(self.p['output_file'])[0], int(self.t))
         with open(restartfile, 'w') as fp:
@@ -2948,7 +2901,6 @@ def parse_callback(self, callback):
         -------
         function
             Python callback function
-
         '''
 
         if isinstance(callback, str):
@@ -2967,18 +2919,17 @@ def parse_callback(self, callback):
         return None
 
 
-    def print_progress(self, fraction=.01, min_interval=1., max_interval=60.):
+    def print_progress(self, fraction:float=.01, min_interval:float=1., max_interval:float=60.) -> None:
         '''Print progress to screen
 
         Parameters
         ----------
-        fraction : float, optional
+        fraction : optional
             Fraction of simulation at which to print progress (default: 1%)
-        min_interval : float, optional
+        min_interval :  optional
             Minimum time in seconds between subsequent progress prints (default: 1s)
-        max_interval : float, optional
+        max_interval : optional
             Maximum time in seconds between subsequent progress prints (default: 60s)
-
         '''
 
         p = (self.t-self.p['tstart']) / (self.p['tstop']-self.p['tstart'])
@@ -3005,7 +2956,7 @@ def print_progress(self, fraction=.01, min_interval=1., max_interval=60.):
             
 
 
-    def print_params(self):
+    def print_params(self) -> None:
         '''Print model configuration parameters to screen'''
 
         maxl = np.max([len(par) for par in self.p.keys()])
@@ -3033,7 +2984,7 @@ def print_params(self):
         logger.info('')
 
 
-    def print_stats(self):
+    def print_stats(self) -> None:
         '''Print model run statistics to screen'''
 
         n_time = self.get_count('time')
@@ -3085,14 +3036,14 @@ class WindGenerator():
     # source:
     # http://www.lutralutra.co.uk/2012/07/02/simulating-a-wind-speed-time-series-in-python/
 
-
+    
     def __init__(self,
-                 mean_speed=9.0,
-                 max_speed=30.0,
-                 dt=60.,
-                 n_states=30,
-                 shape=2.,
-                 scale=2.):
+                 mean_speed:float=9.0,
+                 max_speed:float=30.0,
+                 dt:float=60., # size of time step
+                 n_states:int=30,
+                 shape:float=2.,
+                 scale:float=2.) -> None:
 
         self.mean_speed=mean_speed
         self.max_speed=max_speed
@@ -3146,11 +3097,26 @@ def __init__(self,
         self.MTMcum = np.cumsum(MTM,1)
 
 
-    def __getitem__(self, s):
+    def __getitem__(self, s:int)-> ndarray:
+        """
+        Retrieves item from list of wind speeds as an array 
+        
+        Parameters:
+        -----------
+        s: 
+            index in a list
+        """
+
         return np.asarray(self.wind_speeds[s])
 
 
-    def generate(self, duration=3600.):
+    def generate(self, duration:float=3600.):
+        """Generates wind velocities
+
+        Parameters:
+        -----------
+            duration: duration of time series in seconds
+        """
 
         # initialise series
         self.state = 0
@@ -3168,7 +3134,11 @@ def generate(self, duration=3600.):
         return self
 
 
-    def update(self):
+    def update(self) -> None:
+        """
+        Generates wind speeds base on ramdom uniform distribution 
+        """
+
         r1 = np.random.uniform(0,1)
         r2 = np.random.uniform(0,1)
 
@@ -3184,14 +3154,24 @@ def update(self):
         self.t += self.dt
 
 
-    def get_time_series(self):
+    def get_time_series(self) -> Any:
+        """
+        returns time series
+        """
         u = np.asarray(self.wind_speeds)
         t = np.arange(len(u)) * self.dt
 
         return t, u
 
 
-    def write_time_series(self, fname):
+    def write_time_series(self, fname:str) -> None:
+        """
+        Saves time series to file
+
+        Parameters:
+        -----------
+            fname: file name to write time series to
+        """
         t, u = self.get_time_series()
         M = np.concatenate((np.asmatrix(t),
                             np.asmatrix(u),
@@ -3201,6 +3181,12 @@ def write_time_series(self, fname):
 
 
     def plot(self):
+        """
+        Creates plot of wind speed vs. time
+        
+        Returns:
+
+        """
         t, u = self.get_time_series()
 
         fig, axs = plt.subplots(figsize=(10,4))
@@ -3216,6 +3202,12 @@ def plot(self):
 
 
     def hist(self):
+        """
+        Creates histogram for wind speeds
+
+        Returns:
+
+        """
         fig, axs = plt.subplots(figsize=(10,4))
         axs.hist(self.wind_speeds, bins=self.bins, normed=True, color='k')
         axs.set_xlabel('wind speed [m/s]')
@@ -3227,6 +3219,7 @@ def hist(self):
 
     @staticmethod
     def weibullpdf(data, scale, shape):
+        # TODO: add description
         return [(shape/scale)
                 * ((x/scale)**(shape-1))
                 * np.exp(-1*(x/scale)**shape)
@@ -3235,4 +3228,5 @@ def weibullpdf(data, scale, shape):
 
     @staticmethod
     def matmult4(m, v):
+        # TODO: add description
         return [reduce(operator.add, map(operator.mul,r,v)) for r in m]
diff --git a/aeolis/tests/test_model.py b/aeolis/tests/test_model.py
new file mode 100644
index 00000000..f5aacfab
--- /dev/null
+++ b/aeolis/tests/test_model.py
@@ -0,0 +1,94 @@
+"""
+This file is part of AeoLiS test suit. It uses the Pytest framework.
+
+To run all tests, use:
+    pytest
+
+To run only the test in this file, use:
+    pytest aeolis/tests/test_utils.py
+
+To run all test in a group (class), use:
+    pytest aeolis/tests/test_utils.py::<TestClass>
+
+To run an specific test, use:
+    - If in a group:
+        pytest aeolis/tests/test_utils.py::<TestClass>::<test-method>
+    
+    - If not in any group:
+        pytest aeolis/tests/test_utils.py::<test-function>
+"""
+
+import logging
+import pytest
+import numpy as np
+from aeolis.model import (StreamFormatter, 
+                          ModelState
+                        )
+
+
+class TestStreamFormatter:
+    """Test the stream formatter for console output"""
+
+    def test_stream_formatter_parent(self):
+        """Test if the stream formatter inherits from the logging.Formatter class"""
+        stream_formatter = StreamFormatter()
+        assert isinstance(stream_formatter, logging.Formatter)
+    
+    def test_stream_formatter_info_level(self):
+        """Test if the stream formatter returns log message for the info level"""
+        logger = logging.getLogger("Test")
+        info_record = logger.makeRecord("Test", logging.INFO, 
+                                        "Test", 1, "This is a message for INFO level",
+                                        None, None)
+        stream_formatter = StreamFormatter()
+        info_message = stream_formatter.format(info_record)
+        
+        warning_record = logger.makeRecord("Test warning", logging.WARNING,
+                                        "Test", 2, "This is a message for WARNING level",
+                                        None, None)
+
+        warning_message = stream_formatter.format(warning_record)
+
+        assert info_message == "This is a message for INFO level"
+        assert warning_message == "WARNING: This is a message for WARNING level"
+
+
+class TestModelState:
+    """Test the model state class"""
+
+    def test_model_initialization(self):
+        """Test if the model state class is initialized properly"""
+        state = ModelState(args=None, kwargs=None)
+
+        assert isinstance(state, ModelState)
+
+    def test_set_variable_and_value(self):
+        """Test if the model state class can set a variable and value, and 
+        it is added to the set of mutable variables
+        """
+        state = ModelState(args=None, kwargs=None)
+        state.__setitem__("variable1", 1)
+
+        assert state["variable1"] == 1
+        assert "variable1" in state.ismutable
+
+    def test_variable_is_set_as_mutable(self):
+        """Test if a variable in the model stated is added to
+          the mutable set
+          """
+        state = ModelState(args=None, kwargs=None)
+        state.__setitem__("variable1", 2)
+
+        state.set_mutable("variable2")
+        assert "variable2" in state.ismutable
+    
+    def test_variable_is_set_as_immutable(self):
+        """Test if a variable in the model stated is removed from
+        the mutable set
+          """
+        state = ModelState(args=None, kwargs=None)
+        state.__setitem__("variable1", 3)
+
+        state.set_immutable("variable2")
+        assert "variable2" not in state.ismutable
+
diff --git a/aeolis/tests/test_utils.py b/aeolis/tests/test_utils.py
index 886fc3cb..86666e4b 100644
--- a/aeolis/tests/test_utils.py
+++ b/aeolis/tests/test_utils.py
@@ -18,9 +18,11 @@
         pytest aeolis/tests/test_utils.py::<test-function>
 """
 
+
+
 from aeolis import utils
 import numpy as np
-
+import pytest
 
 # Use test classes to group test-cases that belong to the 
 # same scenario; unless a single test-case is required
@@ -29,6 +31,7 @@ class TestIsIterable:
     """
     Test if iterable variables can be correctly distinguished from non-iterable variables
     """
+    
     def test_isiterable(self):
         x_iterable = np.random.rand(3,2)
         """
@@ -41,4 +44,71 @@ def test_not_iterable(self):
         x_iterable = None
         assert utils.isiterable(x_iterable) == False
 
+class TestMakeIterable:
+
+    """Test if no-iterable variables are correctly returned as iterable variables"""
+
+    def test_make_iterable(self):
+        """Check if a non-iterable object as input returns an iterable object"""
+        x_iterable = None
+        assert utils.isiterable(x_iterable) == False
+        assert utils.isiterable(utils.makeiterable(x_iterable)) == True
+
+class TestIsArray:
+    """ Test if array variables can be correctly distinguished from non-array variables"""
+
+    def test_is_array(self):
+        """Check if an array object as input returns True"""
+        x_array = np.random.rand(3,2)
+        assert utils.isarray(x_array) == True
+
+    def test_not_array(self):
+        """Check if non-array object as input returns False"""
+        x_array = "Hello World!"
+        assert utils.isarray(x_array) == False
+
+@pytest.fixture
+def y():
+    """Return a 2D array"""
+    return np.random.rand(3,3)
+
+class TestInterpArray:
+    """Test if interpolation of multiple time series is correctly performed"""
+
+    def test_interp_array(self, y):
+        pass
+        #TODO: Fix this test
+        # """Check if interpolation of multiple time series is correctly performed"""
+        # x = np.random.rand(3,2)
+        # y = np.random.rand(3,2)
+        
+        # assert utils.isarray(utils.interp_array(x, y, x_new)) == True
+
+
+class TestInterpCircular:
+    """Test the interp_circular function"""
+
+    def test_value_error_raised(self):
+        """Test if a value error is raised when xp and f have different lengths"""
+        
+        x = np.random.rand(10)
+        xp = np.random.rand(10)
+        fp = np.random.rand(9)        
+
+        with pytest.raises(ValueError):
+            utils.interp_circular(x, xp, fp)
+
+
+class TestPreventTinyNegatives:
+    """Test the prevent tiny negatives function"""
+
+    def test_prevent_tiny_negatives(self):
+        """Test if tiny negative values in an array are replace a the replacement value"""
+        
+        x = np.array([1, 2, -1e-10, 3, -1e-10, 4, 5])
+        x = utils.prevent_tiny_negatives(x, max_error=1e-9, replacement=-0.01)
+
+        assert np.all(x >= -0.01)
+
+
 # TODO: Implement other unit-tests.
diff --git a/aeolis/utils.py b/aeolis/utils.py
index 4b38daf9..d886111c 100644
--- a/aeolis/utils.py
+++ b/aeolis/utils.py
@@ -30,6 +30,7 @@
 #import numba
 #import numba_scipy
 import scipy
+from numpy import ndarray
 
 
 def isiterable(x):
@@ -66,7 +67,8 @@ def isarray(x):
         return False
 
 
-def interp_array(x, xp, fp, circular=False, **kwargs):
+def interp_array(x: ndarray, xp: ndarray, 
+                 fp: ndarray, circular: bool=False, **kwargs: dict) -> ndarray:
     '''Interpolate multiple time series at once
 
     Parameters
@@ -99,7 +101,7 @@ def interp_array(x, xp, fp, circular=False, **kwargs):
     return f
 
 
-def interp_circular(x, xp, fp, **kwargs):
+def interp_circular(x: ndarray, xp: ndarray, fp: ndarray, **kwargs) -> ndarray:
     '''One-dimensional linear interpolation.
 
     Returns the one-dimensional piecewise linear interpolant to a
@@ -139,7 +141,7 @@ def interp_circular(x, xp, fp, **kwargs):
     return np.interp(x, xp, fp, **kwargs)
 
 
-def normalize(x, ref=None, axis=0, fill=0.):
+def normalize(x: ndarray, ref:ndarray = None, axis: int =0, fill: float =0.):
     '''Normalize array
 
     Normalizes an array to make it sum to unity over a specific
@@ -152,9 +154,9 @@ def normalize(x, ref=None, axis=0, fill=0.):
         The array to be normalized
     ref : array_like, optional
         Alternative normalization reference, if not specified, the sum of x is used
-    axis : int, optional
+    axis : optional
         The normalization axis (default: 0)
-    fill : float, optional
+    fill : optional
         The return value for all-zero dimensions (default: 0.)
 
     '''
@@ -168,7 +170,7 @@ def normalize(x, ref=None, axis=0, fill=0.):
     return y
 
 
-def prevent_tiny_negatives(x, max_error=1e-10, replacement=0.):
+def prevent_tiny_negatives(x: ndarray, max_error: float =1e-10, replacement: float =0.) -> ndarray:
     '''Replace tiny negative values in array
     
     Parameters