diff --git a/example_PNP.py b/example_PNP.py
index 392f0dc6a8..82cc8b29cb 100644
--- a/example_PNP.py
+++ b/example_PNP.py
@@ -1,9 +1,8 @@
-"""
 import openpnm as op
 import numpy as np
 ws = op.Workspace()
 proj = ws.new_project()
-# ws.settings['loglevel'] = 20
+
 
 # network, geometry, phase
 np.random.seed(0)
@@ -14,10 +13,6 @@
        net['pore.front'] * net['pore.left'])
 thrts = net['throat.surface']
 op.topotools.trim(network=net, pores=net.Ps[prs], throats=net.Ts[thrts])
-net['pore.coords'][:, 0] -= 0.00045
-net['pore.coords'][:, 1] -= 0.00045
-net['pore.coords'][:, 2] = 0
-# np.savetxt('coords.txt', net['pore.coords'])
 
 
 geo = op.geometry.StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
@@ -28,8 +23,6 @@
 geo.regenerate_models()
 
 sw = op.phases.Mixtures.SalineWater(network=net)
-sw['pore.diffusivity.Na'] = sw['pore.diffusivity.Cl']
-sw['throat.diffusivity.Na'] = sw['throat.diffusivity.Cl']
 
 # physics
 phys = op.physics.GenericPhysics(network=net, phase=sw, geometry=geo)
@@ -67,17 +60,17 @@
 p = op.algorithms.OhmicConduction(network=net, phase=sw)
 p.settings['conductance'] = 'throat.ionic_conductance'
 p.settings['quantity'] = 'pore.potential'
-p.set_value_BC(pores=net.pores('left'), values=0.05)
+p.set_value_BC(pores=net.pores('left'), values=0.01)
 p.set_value_BC(pores=net.pores('right'), values=0.00)
 p.settings['rxn_tolerance'] = 1e-12
-sw.update(p.results())
 
-eA = op.algorithms.NernstPlanck(network=net, phase=sw, electrolyte='Na')
+
+eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion='Na')
 eA.set_value_BC(pores=net.pores('back'), values=100)
 eA.set_value_BC(pores=net.pores('front'), values=90)
 eA.settings['rxn_tolerance'] = 1e-12
 
-eB = op.algorithms.NernstPlanck(network=net, phase=sw, electrolyte='Cl')
+eB = op.algorithms.NernstPlanck(network=net, phase=sw, ion='Cl')
 eB.set_value_BC(pores=net.pores('back'), values=100)
 eB.set_value_BC(pores=net.pores('front'), values=90)
 eB.settings['rxn_tolerance'] = 1e-12
@@ -93,27 +86,14 @@
                model=ad_dif_mig_Cl, ion='Cl',
                s_scheme='exponential')
 
-pnp = op.algorithms.PoissonNernstPlanck(network=net, phase=sw)
-pnp.setup(potential_field=p, electrolytes=[eA, eB])
-pnp.settings['max_iter'] = 400
+pnp = op.algorithms.ChargeConservationNernstPlanck(network=net, phase=sw)
+pnp.setup(potential_field=p, ions=[eA, eB])
+pnp.settings['max_iter'] = 10
 pnp.settings['tolerance'] = 1e-04
-# pnp.settings['charge_conservation'] = 'laplace'
-pnp.settings['charge_conservation'] = 'electroneutrality'
+pnp.settings['charge_conservation'] = 'laplace'
 pnp.run()
 
-# Comsol results
-# cNa_cmsl = np.genfromtxt('cNa_at_coords.txt', skip_header=9)
-# sw['pore.concentration.cmsl'] = cNa_cmsl[:, 3]
-
-# phi_cmsl = np.genfromtxt('phi_at_coords.txt', skip_header=9)
-# sw['pore.potential_cmsl'] = phi_cmsl[:, 3]
-
-# p_cmsl = np.genfromtxt('p_at_coords.txt', skip_header=9)
-# sw['pore.pressure_cmsl'] = p_cmsl[:, 3]
-
 sw.update(sf.results())
 sw.update(p.results())
 sw.update(eA.results())
 sw.update(eB.results())
-# op.io.VTK.save(network=net, phases=[sw], filename='OUTPUT_PNP')
-"""
diff --git a/example_ad_01.py b/example_ad_01.py
deleted file mode 100644
index 53ba4295a8..0000000000
--- a/example_ad_01.py
+++ /dev/null
@@ -1,56 +0,0 @@
-import openpnm as op
-import numpy as np
-
-# work space and project
-ws = op.Workspace()
-proj = ws.new_project()
-
-# network
-np.random.seed(7)
-n = 10
-net = op.network.Cubic(shape=[n, n, 1], spacing=1e-4, project=proj)
-
-# geometry
-geo = op.geometry.StickAndBall(network=net,
-                               pores=net.Ps,
-                               throats=net.Ts)
-
-# phase
-phase = op.phases.Water(network=net)
-
-# physics
-phys = op.physics.GenericPhysics(network=net,
-                                 phase=phase,
-                                 geometry=geo)
-phase['pore.diffusivity'] = 2e-09
-phase['throat.diffusivity'] = 2e-09
-
-mod1 = op.models.physics.hydraulic_conductance.hagen_poiseuille
-phys.add_model(propname='throat.hydraulic_conductance',
-               model=mod1, regen_mode='normal')
-
-# algorithms: Stokes flow
-sf = op.algorithms.StokesFlow(network=net, phase=phase)
-sf.set_value_BC(pores=net.pores('left'), values=0)
-sf.set_value_BC(pores=net.pores('right'), values=0.1)
-sf.run()
-phase.update(sf.results())
-
-
-# algorithms: dispersion
-mod2 = op.models.physics.diffusive_conductance.ordinary_diffusion
-phys.add_model(propname='throat.diffusive_conductance',
-               model=mod2, regen_mode='normal')
-
-mod3 = op.models.physics.dispersive_conductance.dispersion
-phys.add_model(propname='throat.dispersive_conductance',
-               model=mod3, regen_mode='normal')
-
-dis = op.algorithms.Dispersion(network=net, phase=phase)
-dis.set_value_BC(pores=net.pores('back'), values=1)
-dis.set_value_BC(pores=net.pores('front'), values=0.5)
-dis.run()
-
-# output results to a vtk file
-# phase.update(dis.results())
-# proj.export_data(phases=[phase], filename='OUT', filetype='VTK')
diff --git a/example_ad_02.py b/example_ad_02.py
deleted file mode 100644
index 499caabee6..0000000000
--- a/example_ad_02.py
+++ /dev/null
@@ -1,63 +0,0 @@
-import openpnm as op
-import numpy as np
-
-# work space and project
-ws = op.Workspace()
-proj = ws.new_project()
-
-# network
-np.random.seed(7)
-n = 10
-net = op.network.Cubic(shape=[n, n, 1], spacing=1e-4, project=proj)
-
-# geometry
-geo = op.geometry.StickAndBall(network=net,
-                               pores=net.Ps,
-                               throats=net.Ts)
-
-# phase
-phase = op.phases.Water(network=net)
-
-# physics
-phys = op.physics.GenericPhysics(network=net,
-                                 phase=phase,
-                                 geometry=geo)
-phase['pore.diffusivity'] = 2e-09
-phase['throat.diffusivity'] = 2e-09
-
-mod1 = op.models.physics.hydraulic_conductance.hagen_poiseuille
-phys.add_model(propname='throat.hydraulic_conductance',
-               model=mod1, regen_mode='normal')
-
-# algorithms: Stokes flow
-sf = op.algorithms.StokesFlow(network=net, phase=phase)
-sf.set_value_BC(pores=net.pores('left'), values=0)
-sf.set_value_BC(pores=net.pores('right'), values=0.1)
-sf.run()
-phase.update(sf.results())
-
-
-# algorithms: dispersion
-mod2 = op.models.physics.diffusive_conductance.ordinary_diffusion
-phys.add_model(propname='throat.diffusive_conductance',
-               model=mod2, regen_mode='normal')
-
-mod3 = op.models.physics.ad_dif_conductance.ad_dif
-phys.add_model(propname='throat.ad_dif_conductance',
-               model=mod3, regen_mode='normal')
-
-dis = op.algorithms.AdvectionDiffusion(network=net, phase=phase)
-dis.set_value_BC(pores=net.pores('back'), values=1)
-dis.set_value_BC(pores=net.pores('front'), values=0.5)
-dis.run()
-
-dis2 = op.algorithms.TransientAdvectionDiffusion(network=net, phase=phase)
-dis2.set_value_BC(pores=net.pores('back'), values=1)
-dis2.set_value_BC(pores=net.pores('front'), values=0.5)
-dis2.setup(t_initial=0, t_final=300, t_step=1, t_output=50, t_tolerance=1e-20,
-           t_scheme='implicit')
-dis2.run()
-
-# output results to a vtk file
-# phase.update(dis.results())
-# proj.export_data(phases=[phase], filename='OUT', filetype='VTK')
diff --git a/openpnm/algorithms/PoissonNernstPlanck.py b/openpnm/algorithms/ChargeConservationNernstPlanck.py
similarity index 77%
rename from openpnm/algorithms/PoissonNernstPlanck.py
rename to openpnm/algorithms/ChargeConservationNernstPlanck.py
index 6e30d66ef1..faade5077f 100644
--- a/openpnm/algorithms/PoissonNernstPlanck.py
+++ b/openpnm/algorithms/ChargeConservationNernstPlanck.py
@@ -5,14 +5,15 @@
 logger = logging.getLogger(__name__)
 
 
-class PoissonNernstPlanck(ReactiveTransport):
+class ChargeConservationNernstPlanck(ReactiveTransport):
     r"""
-    A subclass of GenericTransport to solve the Poisson Nernst-Planck equations
+    A subclass of GenericTransport to solve the charge conservation and
+    Nernst-Planck equations.
     """
     def __init__(self, settings={}, phase=None, **kwargs):
         def_set = {'phase': None,
                    'potential_field': None,
-                   'electrolytes': [],
+                   'ions': [],
                    'charge_conservation': 'electroneutrality',
                    'tolerance': 1e-4,
                    'max_iter': 10}
@@ -22,7 +23,7 @@ def __init__(self, settings={}, phase=None, **kwargs):
         if phase is not None:
             self.setup(phase=phase)
 
-    def setup(self, phase=None, potential_field=None, electrolytes=[],
+    def setup(self, phase=None, potential_field=None, ions=[],
               charge_conservation=None, tolerance=None,
               max_iter=None, **kwargs):
         r"""
@@ -31,8 +32,8 @@ def setup(self, phase=None, potential_field=None, electrolytes=[],
             self.settings['phase'] = phase.name
         if potential_field:
             self.settings['potential_field'] = potential_field
-        if electrolytes:
-            self.settings['electrolytes'] = electrolytes
+        if ions:
+            self.settings['ions'] = ions
         if charge_conservation:
             self.settings['charge_conservation'] = charge_conservation
         if tolerance:
@@ -42,25 +43,32 @@ def setup(self, phase=None, potential_field=None, electrolytes=[],
         super().setup(**kwargs)
 
     def run(self):
-        # Phase, potential and electrolytes algorithms
+        # Phase, potential and ions algorithms
         phase = self.project.phases()[self.settings['phase']]
         p_alg = self.settings['potential_field']
-        e_alg = self.settings['electrolytes']
+        e_alg = self.settings['ions']
 
         # Define initial conditions (if not defined by the user)
         try:
             p_alg[p_alg.settings['quantity']]
         except KeyError:
-            p_alg[p_alg.settings['quantity']] = np.zeros(shape=[p_alg.Np, ],
-                                                         dtype=float)
+            try:
+                p_alg[p_alg.settings['quantity']] = (
+                    phase[p_alg.settings['quantity']])
+            except KeyError:
+                p_alg[p_alg.settings['quantity']] = np.zeros(
+                    shape=[p_alg.Np, ], dtype=float)
         for e in e_alg:
             try:
                 e[e.settings['quantity']]
             except KeyError:
-                e[e.settings['quantity']] = np.zeros(shape=[e.Np, ],
-                                                     dtype=float)
+                try:
+                    e[e.settings['quantity']] = phase[e.settings['quantity']]
+                except KeyError:
+                    e[e.settings['quantity']] = np.zeros(shape=[e.Np, ],
+                                                         dtype=float)
 
-        # Define source term for Poisson equation
+        # Source term for Poisson or charge conservation (electroneutrality) eq
         Ps = (p_alg['pore.all'] * np.isnan(p_alg['pore.bc_value']) *
               np.isnan(p_alg['pore.bc_rate']))
         mod = gst.charge_conservation
@@ -80,8 +88,7 @@ def run(self):
         # Iterate until solutions converge
         for itr in range(int(self.settings['max_iter'])):
             r = str([float(format(i, '.3g')) for i in res.values()])[1:-1]
-            logger.info('Iter: ' + str(itr) + ', Residuals: ' + r)
-            print('Iter: ' + str(itr) + ', Residuals: ' + r)
+            logger.info('Iter: ' + str(itr+1) + ', Residuals: ' + r)
             convergence = max(i for i in res.values()) < tol
             if not convergence:
                 # Poisson eq
@@ -95,7 +102,7 @@ def run(self):
                 phase.update(p_alg.results())
                 phys[0].regenerate_models()
 
-                # Electrolytes
+                # Ions
                 for e in e_alg:
                     c_old = e[e.settings['quantity']].copy()
                     e._run_reactive(x=c_old)
diff --git a/openpnm/algorithms/NernstPlanck.py b/openpnm/algorithms/NernstPlanck.py
index cc4509625c..01a44ce303 100644
--- a/openpnm/algorithms/NernstPlanck.py
+++ b/openpnm/algorithms/NernstPlanck.py
@@ -5,13 +5,16 @@
 
 class NernstPlanck(ReactiveTransport):
     r"""
+    A class to simulate transport of charged species (such as ions) in dilute
+    solutions.
+
     """
-    def __init__(self, settings={}, phase=None, electrolyte='', **kwargs):
+    def __init__(self, settings={}, phase=None, ion='', **kwargs):
         def_set = {'phase': None,
-                   'quantity': 'pore.concentration.'+electrolyte,
-                   'conductance': 'throat.ad_dif_mig_conductance.'+electrolyte}
+                   'quantity': 'pore.concentration.'+ion,
+                   'conductance': 'throat.ad_dif_mig_conductance.'+ion}
         super().__init__(**kwargs)
-        self.name = electrolyte
+        self.name = ion
         self.settings.update(def_set)
         self.settings.update(settings)
         if phase is not None:
diff --git a/openpnm/algorithms/__init__.py b/openpnm/algorithms/__init__.py
index cd7c566d24..81163fe3a1 100644
--- a/openpnm/algorithms/__init__.py
+++ b/openpnm/algorithms/__init__.py
@@ -27,4 +27,4 @@
 from .MixedInvasionPercolation import MixedInvasionPercolation
 from .Porosimetry import Porosimetry
 from .NernstPlanck import NernstPlanck
-from .PoissonNernstPlanck import PoissonNernstPlanck
+from .ChargeConservationNernstPlanck import ChargeConservationNernstPlanck
diff --git a/openpnm/models/physics/ad_dif_mig_conductance.py b/openpnm/models/physics/ad_dif_mig_conductance.py
index 41848f7e8a..170cd0aa14 100644
--- a/openpnm/models/physics/ad_dif_mig_conductance.py
+++ b/openpnm/models/physics/ad_dif_mig_conductance.py
@@ -233,9 +233,15 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
             T = phase[throat_temperature][throats]
         except KeyError:
             T = phase.interpolate_data(propname=pore_temperature)[throats]
-
-        P = phase[pore_pressure]
-        V = phase[pore_potential]
+        # Check if pressure and potential values exist, otherwise, assign zeros
+        try:
+            P = phase[pore_pressure]
+        except KeyError:
+            P = _sp.zeros(shape=[phase.Np, ], dtype=float)
+        try:
+            V = phase[pore_potential]
+        except KeyError:
+            V = _sp.zeros(shape=[phase.Np, ], dtype=float)
         gh = phase[throat_hydraulic_conductance]
         gd = phase[throat_diffusive_conductance]
         gd = _sp.tile(gd, 2)
@@ -244,16 +250,19 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
         F = 96485.3329
         R = 8.3145
 
-        S = (A1*L1+A2*L2+At*Lt)/(L1+L2+Lt)
-        L = L1 + Lt + L2
-
         # Advection
         Qij = -gh*_sp.diff(P[cn], axis=1).squeeze()
         Qij = _sp.append(Qij, -Qij)
 
         # Migration
-        grad_V = _sp.diff(V[cn], axis=1).squeeze() / L
-        mig = ((z*F*D*S)/(R*T)) * grad_V
+        gm1, gm2, gmt = _sp.zeros((3, len(Lt)))
+        gm1[~m1] = gm2[~m2] = gmt[~mt] = _sp.inf
+        gm1[m1] = ((z*F*D*A1)/(R*T))[m1] / L1[m1]
+        gm2[m2] = ((z*F*D*A2)/(R*T))[m2] / L2[m2]
+        gmt[mt] = ((z*F*D*At)/(R*T))[mt] / Lt[mt]
+        gm = (1/gm1 + 1/gm2 + 1/gmt)**(-1)
+        delta_V = _sp.diff(V[cn], axis=1).squeeze()
+        mig = gm * delta_V
         mig = _sp.append(mig, -mig)
 
         # Advection-migration
diff --git a/openpnm/models/physics/diffusive_conductance.py b/openpnm/models/physics/diffusive_conductance.py
index 54afb9116f..b4c5a0efbd 100644
--- a/openpnm/models/physics/diffusive_conductance.py
+++ b/openpnm/models/physics/diffusive_conductance.py
@@ -337,7 +337,7 @@ def classic_ordinary_diffusion(target,
     # Get shape factor
     try:
         sf = network[shape_factor]
-    except:
+    except KeyError:
         sf = _sp.ones(network.num_throats())
     sf[_sp.isnan(sf)] = 1.0
     gt = (1 / sf) * ct * DABt * tarea / tlen
diff --git a/openpnm/models/physics/generic_source_term.py b/openpnm/models/physics/generic_source_term.py
index 08dd779e88..bb9ffa4a64 100644
--- a/openpnm/models/physics/generic_source_term.py
+++ b/openpnm/models/physics/generic_source_term.py
@@ -1,5 +1,6 @@
 r"""
 
+.. autofunction:: openpnm.models.physics.generic_source_term.charge_conservation
 .. autofunction:: openpnm.models.physics.generic_source_term.standard_kinetics
 .. autofunction:: openpnm.models.physics.generic_source_term.linear
 .. autofunction:: openpnm.models.physics.generic_source_term.power_law
@@ -17,6 +18,44 @@
 
 
 def charge_conservation(target, phase, p_alg, e_alg, assumption):
+    r"""
+    Applies the source term on the charge conservation equation when solving
+    for ions transport.
+
+    Parameters
+    ----------
+    phase : OpenPNM Phase object
+            The phase on which the charge conservation equation is applied.
+
+    p_alg : OpenPNM Algorithm object
+            The algorithm used to enforce charge conservation.
+
+    e_alg : list of OpenPNM algorithms
+            The list of algorithms used to solve for transport of different
+            ionic species of the mixture phase.
+
+    assumption : string
+            A string correponding to the assumption adopted to enforce charge
+            conservation.
+
+    Returns
+    -------
+    A dictionary containing the following three items:
+
+        **'rate'** - The value of the source term function for the given list
+                     of algortihms under the provided assumption.
+
+        **'S1'** - A placeholder (zero array).
+
+        **'S2'** - The value of the source term function for the given list of
+                   algortihms under the provided assumption (same as 'rate').
+
+    Notes
+    -----
+    Three assumptions are supported; "poisson", "electroneutrality" and
+    "laplace".
+
+    """
     F = 96485.3329
     rhs = _sp.zeros(shape=(p_alg.Np, ), dtype=float)
     if assumption == 'poisson':
diff --git a/openpnm/models/physics/ionic_conductance.py b/openpnm/models/physics/ionic_conductance.py
index b2dba8cf6c..c71a5072df 100644
--- a/openpnm/models/physics/ionic_conductance.py
+++ b/openpnm/models/physics/ionic_conductance.py
@@ -13,8 +13,16 @@
 def poisson(target,
             pore_area='pore.area',
             throat_area='throat.area',
+            pore_diffusivity='pore.diffusivity',
+            throat_diffusivity='throat.diffusivity',
             conduit_lengths='throat.conduit_lengths',
-            conduit_shape_factors='throat.poisson_shape_factors'):
+            conduit_shape_factors='throat.poisson_shape_factors',
+            pore_volume='pore.volume',
+            pore_temperature='pore.temperature',
+            throat_temperature='throat.temperature',
+            pore_valence='pore.valence',
+            throat_valence='throat.valence',
+            pore_concentration='pore.concentration'):
     r"""
     Calculate the ionic conductance of conduits in network (using the Poisson
     equation for charge conservation), where a conduit is
@@ -33,12 +41,36 @@ def poisson(target,
     throat_area : string
         Dictionary key of the throat area values
 
+    pore_diffusivity : string
+        Dictionary key of the pore diffusivity values
+
+    throat_diffusivity : string
+        Dictionary key of the throat diffusivity values
+
     conduit_lengths : string
         Dictionary key of the conduit length values
 
     conduit_shape_factors : string
         Dictionary key of the conduit DIFFUSION shape factor values
 
+    pore_volume : string
+        Dictionary key of the pore volume values
+
+    pore_temperature : string
+        Dictionary key of the pore temperature values
+
+    throat_temperature : string
+        Dictionary key of the throat temperature values
+
+    pore_valence : string
+       Dictionary key of the pore ionic species valence values
+
+    throat_valence : string
+       Dictionary key of the throat ionic species valence values
+
+    pore_concentration : string
+       Dictionary key of the pore ionic species concentration values
+
     Returns
     -------
     g : ndarray
@@ -62,10 +94,15 @@ def poisson(target,
                                transport_type='poisson',
                                pore_area=pore_area,
                                throat_area=throat_area,
-                               pore_diffusivity='',
-                               throat_diffusivity='',
+                               pore_diffusivity=pore_diffusivity,
+                               throat_diffusivity=throat_diffusivity,
                                conduit_lengths=conduit_lengths,
-                               pore_volume='',
+                               pore_volume=pore_volume,
+                               pore_temperature=pore_temperature,
+                               throat_temperature=throat_temperature,
+                               pore_valence=pore_valence,
+                               throat_valence=throat_valence,
+                               pore_concentration=pore_concentration,
                                conduit_shape_factors=conduit_shape_factors)
 
 
@@ -367,9 +404,9 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
         SF1 = SF2 = SFt = 1.0
     # Poisson or Laplace
     if transport_type in ['poisson', 'laplace']:
-        g1[m1] = (A1)[m1] * L1[m1]
-        g2[m2] = (A2)[m2] * L2[m2]
-        gt[mt] = (At)[mt] * Lt[mt]
+        g1[m1] = (A1)[m1] / L1[m1]
+        g2[m2] = (A2)[m2] / L2[m2]
+        gt[mt] = (At)[mt] / Lt[mt]
     # Electroneutrality
     elif transport_type == 'electroneutrality':
         F = 96485.3329
@@ -402,7 +439,7 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
                 c1 = _sp.zeros((network.Nt))[cn[:, 0]]
                 c2 = _sp.zeros((network.Nt))[cn[:, 1]]
             ct = (c1*Vol1 + c2*Vol2)/(Vol1 + Vol2)
-        # Interpolate pore phase property values to throats
+            # Interpolate pore phase property values to throats
             try:
                 Dt = phase[throat_diffusivity+i][throats]
                 Vt = phase[throat_valence+i][throats]
@@ -436,8 +473,11 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
     # Preallocating g_inv
     g_inv1, g_inv2, g_invt = _sp.zeros((3, len(Lt)))
     f1, f2, ft = [gi != 0 for gi in [g1, g2, gt]]
+    g_inv1[~f1] = g_inv2[~f2] = g_invt[~ft] = _sp.inf
     g_inv1[f1] = 1/g1[f1]
     g_inv2[f2] = 1/g2[f2]
     g_invt[ft] = 1/gt[ft]
+    g = g_inv1/SF1 + g_inv2/SF2 + g_invt/SFt
+    g[g != 0] = g[g != 0]**(-1)
     # Apply shape factors and calculate the final conductance
-    return (g_inv1/SF1 + g_inv2/SF2 + g_invt/SFt)**(-1)
+    return g
diff --git a/tests/unit/algorithms/AdvectionDiffusionTest.py b/tests/unit/algorithms/AdvectionDiffusionTest.py
index 186e3b360c..b181671adc 100644
--- a/tests/unit/algorithms/AdvectionDiffusionTest.py
+++ b/tests/unit/algorithms/AdvectionDiffusionTest.py
@@ -34,7 +34,7 @@ def setup_class(self):
         self.ad.set_value_BC(pores=self.net.pores('back'), values=2)
         self.ad.set_value_BC(pores=self.net.pores('front'), values=0)
 
-    def test_powerlaw_advection_diffusion_diffusion(self):
+    def test_powerlaw_advection_diffusion(self):
         mod = op.models.physics.ad_dif_conductance.ad_dif
         self.phys.add_model(propname='throat.ad_dif_conductance_powerlaw',
                             model=mod, s_scheme='powerlaw')
@@ -49,7 +49,7 @@ def test_powerlaw_advection_diffusion_diffusion(self):
         y = sp.around(self.ad['pore.concentration'], decimals=5)
         assert_allclose(actual=y, desired=x)
 
-    def test_upwind_advection_diffusion_diffusion(self):
+    def test_upwind_advection_diffusion(self):
         mod = op.models.physics.ad_dif_conductance.ad_dif
         self.phys.add_model(propname='throat.ad_dif_conductance_upwind',
                             model=mod, s_scheme='upwind')
@@ -63,7 +63,7 @@ def test_upwind_advection_diffusion_diffusion(self):
         y = sp.around(self.ad['pore.concentration'], decimals=5)
         assert_allclose(actual=y, desired=x)
 
-    def test_hybrid_advection_diffusion_diffusion(self):
+    def test_hybrid_advection_diffusion(self):
         mod = op.models.physics.ad_dif_conductance.ad_dif
         self.phys.add_model(propname='throat.ad_dif_conductance_hybrid',
                             model=mod, s_scheme='hybrid')
@@ -78,7 +78,7 @@ def test_hybrid_advection_diffusion_diffusion(self):
         y = sp.around(self.ad['pore.concentration'], decimals=5)
         assert_allclose(actual=y, desired=x)
 
-    def test_exponential_advection_diffusion_diffusion(self):
+    def test_exponential_advection_diffusion(self):
         mod = op.models.physics.ad_dif_conductance.ad_dif
         self.phys.add_model(propname='throat.ad_dif_conductance_exponential',
                             model=mod, s_scheme='exponential')
diff --git a/tests/unit/algorithms/NernstPlanckTest.py b/tests/unit/algorithms/NernstPlanckTest.py
new file mode 100644
index 0000000000..ce6ada502f
--- /dev/null
+++ b/tests/unit/algorithms/NernstPlanckTest.py
@@ -0,0 +1,124 @@
+import openpnm as op
+import scipy as sp
+from numpy.testing import assert_allclose
+
+
+class NernstPlanckTest:
+
+    def setup_class(self):
+        sp.random.seed(0)
+        self.net = op.network.Cubic(shape=[4, 3, 1], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['pore.area'] = 2e-14
+        self.geo['throat.area'] = 1e-14
+        self.geo['throat.conduit_lengths.pore1'] = 0.1
+        self.geo['throat.conduit_lengths.throat'] = 0.6
+        self.geo['throat.conduit_lengths.pore2'] = 0.1
+
+        self.phase = op.phases.GenericPhase(network=self.net)
+        self.phase['pore.diffusivity.ionX'] = 1e-9
+        self.phase['throat.valence.ionX'] = 1
+
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+        self.phys['throat.diffusive_conductance.ionX'] = 1e-15
+        self.phys['throat.hydraulic_conductance'] = 1e-15
+        self.phys['throat.ionic_conductance'] = 1e-15
+
+        self.sf = op.algorithms.StokesFlow(network=self.net, phase=self.phase)
+        self.sf.set_value_BC(pores=self.net.pores('back'), values=1)
+        self.sf.set_value_BC(pores=self.net.pores('front'), values=0)
+        self.sf.run()
+
+        self.phase.update(self.sf.results())
+
+        self.p = op.algorithms.OhmicConduction(network=self.net,
+                                               phase=self.phase)
+        self.p.settings['conductance'] = 'throat.ionic_conductance'
+        self.p.settings['quantity'] = 'pore.potential'
+        self.p.set_value_BC(pores=self.net.pores('left'), values=0.05)
+        self.p.set_value_BC(pores=self.net.pores('right'), values=0.00)
+        self.p.run()
+
+        self.phase.update(self.p.results())
+
+        self.adm = op.algorithms.NernstPlanck(network=self.net,
+                                              phase=self.phase, ion='ionX')
+        self.adm.set_value_BC(pores=self.net.pores('back'), values=2)
+        self.adm.set_value_BC(pores=self.net.pores('front'), values=0)
+
+    def test_powerlaw_NernstPlanck(self):
+        mod = op.models.physics.ad_dif_mig_conductance.ad_dif_mig
+        self.phys.add_model(propname='throat.ad_dif_mig_conductance_powerlaw',
+                            model=mod, s_scheme='powerlaw', ion='ionX')
+        self.phys.regenerate_models()
+
+        self.adm.setup(conductance='throat.ad_dif_mig_conductance_powerlaw')
+        self.adm.run()
+        x = [0.,      0.,      0.,
+             0.89653, 0.89653, 0.89653,
+             1.53924, 1.53924, 1.53924,
+             2.,      2.,      2.]
+        y = sp.around(self.adm['pore.concentration.ionX'], decimals=5)
+        assert_allclose(actual=y, desired=x)
+
+    def test_upwind_NernstPlanck(self):
+        mod = op.models.physics.ad_dif_mig_conductance.ad_dif_mig
+        self.phys.add_model(propname='throat.ad_dif_mig_conductance_upwind',
+                            model=mod, s_scheme='upwind', ion='ionX')
+        self.phys.regenerate_models()
+        self.adm.setup(conductance='throat.ad_dif_mig_conductance_upwind')
+        self.adm.run()
+        x = [0.,      0.,      0.,
+             0.86486, 0.86486, 0.86486,
+             1.51351, 1.51351, 1.51351,
+             2.,      2.,      2.]
+        y = sp.around(self.adm['pore.concentration.ionX'], decimals=5)
+        assert_allclose(actual=y, desired=x)
+
+    def test_hybrid_NernstPlanck(self):
+        mod = op.models.physics.ad_dif_mig_conductance.ad_dif_mig
+        self.phys.add_model(propname='throat.ad_dif_mig_conductance_hybrid',
+                            model=mod, s_scheme='hybrid', ion='ionX')
+        self.phys.regenerate_models()
+
+        self.adm.setup(conductance='throat.ad_dif_mig_conductance_hybrid')
+        self.adm.run()
+        x = [0.,      0.,      0.,
+             0.89908, 0.89908, 0.89908,
+             1.54128, 1.54128, 1.54128,
+             2.,      2.,      2.]
+        y = sp.around(self.adm['pore.concentration.ionX'], decimals=5)
+        assert_allclose(actual=y, desired=x)
+
+    def test_exponential_NernstPlanck(self):
+        mod = op.models.physics.ad_dif_mig_conductance.ad_dif_mig
+        self.phys.add_model(propname='throat.ad_dif_mig_conductance_exp',
+                            model=mod, s_scheme='exponential', ion='ionX')
+        self.phys.regenerate_models()
+
+        self.adm.setup(conductance='throat.ad_dif_mig_conductance_exp')
+        self.adm.run()
+        x = [0.,      0.,      0.,
+             0.89688, 0.89688, 0.89688,
+             1.53953, 1.53953, 1.53953,
+             2.,      2.,      2.]
+        y = sp.around(self.adm['pore.concentration.ionX'], decimals=5)
+        assert_allclose(actual=y, desired=x)
+
+    def teardown_class(self):
+        ws = op.Workspace()
+        ws.clear()
+
+
+if __name__ == '__main__':
+    t = NernstPlanckTest()
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
+    self = t