Repository: PMEAL/OpenPNM
Pull request #1272 (2019-05-05 14:04:54)

Enhancing mixtures
The new Mixture class was not well thought out, so this PR attempts to 'complete' it.  Things have been re-organized in the subfolders, and new methods have been added to the GenericMixture class, among other changes.

Note that this PR has @magnaou's PNP algorithm merged into it, so I could ensure that my SalineWater worked with his example script.  


Patch:
diff --git a/example_PNP.py b/example_PNP.py
index 82cc8b29c..1efb4ba3b 100644
--- a/example_PNP.py
+++ b/example_PNP.py
@@ -1,7 +1,9 @@
 import openpnm as op
+from openpnm.phases import mixtures
 import numpy as np
 ws = op.Workspace()
 proj = ws.new_project()
+# ws.settings['loglevel'] = 20
 
 
 # network, geometry, phase
@@ -22,7 +24,9 @@ geo.add_model(propname='pore.diameter', model=pore_d, value=1.5e-4)
 geo.add_model(propname='throat.diameter', model=throat_d, value=1e-4)
 geo.regenerate_models()
 
-sw = op.phases.Mixtures.SalineWater(network=net)
+sw = mixtures.SalineWater(network=net)
+# Retrieve handles to each species for use below
+Cl, Na, H2O = sw.components.values()
 
 # physics
 phys = op.physics.GenericPhysics(network=net, phase=sw, geometry=geo)
@@ -38,15 +42,15 @@ phys.add_model(propname='throat.ionic_conductance',
                model=current, regen_mode='normal')
 
 eA_dif = op.models.physics.diffusive_conductance.ordinary_diffusion
-phys.add_model(propname='throat.diffusive_conductance.Na',
-               pore_diffusivity='pore.diffusivity.Na',
-               throat_diffusivity='throat.diffusivity.Na',
+phys.add_model(propname='throat.diffusive_conductance.' + Na.name,
+               pore_diffusivity='pore.diffusivity.' + Na.name,
+               throat_diffusivity='throat.diffusivity.' + Na.name,
                model=eA_dif, regen_mode='normal')
 
 eB_dif = op.models.physics.diffusive_conductance.ordinary_diffusion
-phys.add_model(propname='throat.diffusive_conductance.Cl',
-               pore_diffusivity='pore.diffusivity.Cl',
-               throat_diffusivity='throat.diffusivity.Cl',
+phys.add_model(propname='throat.diffusive_conductance.' + Cl.name,
+               pore_diffusivity='pore.diffusivity.' + Cl.name,
+               throat_diffusivity='throat.diffusivity.' + Cl.name,
                model=eB_dif, regen_mode='normal')
 
 # algorithms
@@ -63,27 +67,28 @@ p.settings['quantity'] = 'pore.potential'
 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
+p.run()
+sw.update(p.results())
 
-
-eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion='Na')
+eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion=Na.name)
 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, ion='Cl')
+eB = op.algorithms.NernstPlanck(network=net, phase=sw, ion=Cl.name)
 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
 
 ad_dif_mig_Na = op.models.physics.ad_dif_mig_conductance.ad_dif_mig
-phys.add_model(propname='throat.ad_dif_mig_conductance.Na',
-               model=ad_dif_mig_Na, ion='Na',
+phys.add_model(propname='throat.ad_dif_mig_conductance.' + Na.name,
+               model=ad_dif_mig_Na, ion=Na.name,
                s_scheme='exponential')
 
 ad_dif_mig_Cl = op.models.physics.ad_dif_mig_conductance.ad_dif_mig
-phys.add_model(propname='throat.ad_dif_mig_conductance.Cl',
+phys.add_model(propname='throat.ad_dif_mig_conductance.' + Cl.name,
                pore_pressure=sf.settings['quantity'],
-               model=ad_dif_mig_Cl, ion='Cl',
+               model=ad_dif_mig_Cl, ion=Cl.name,
                s_scheme='exponential')
 
 pnp = op.algorithms.ChargeConservationNernstPlanck(network=net, phase=sw)
@@ -91,6 +96,8 @@ pnp.setup(potential_field=p, ions=[eA, eB])
 pnp.settings['max_iter'] = 10
 pnp.settings['tolerance'] = 1e-04
 pnp.settings['charge_conservation'] = 'laplace'
+# Electroneutrality condition does not work with new Mixtures
+# pnp.settings['charge_conservation'] = 'electroneutrality'
 pnp.run()
 
 sw.update(sf.results())
diff --git a/example_script_mixtures.py b/example_script_mixtures.py
new file mode 100644
index 000000000..abece19da
--- /dev/null
+++ b/example_script_mixtures.py
@@ -0,0 +1,48 @@
+import openpnm as op
+from openpnm.phases import mixtures
+ws = op.Workspace()
+proj = ws.new_project()
+
+pn = op.network.Cubic(shape=[30, 30, 10], spacing=1e-4, project=proj)
+geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
+
+N2 = mixtures.species.gases.N2(network=pn, name='pure_N2')
+O2 = mixtures.species.gases.O2(network=pn, name='pure_O2')
+air = mixtures.GenericMixture(network=pn, components=[N2, O2],
+                              name='air_mixture')
+air.set_mole_fraction(N2, 0.79)
+air.set_mole_fraction(O2, 0.21)
+air.add_model(propname='pore.molar_mass',
+              model=op.models.phases.mixtures.mole_weighted_average,
+              prop='pore.molecular_weight')
+air.add_model(propname='pore.diffusivity.pure_O2',
+              model=op.models.phases.mixtures.fuller_diffusivity)
+air.add_model(propname='pore.viscosity',
+              model=op.models.misc.polynomial,
+              prop='pore.temperature',
+              a=[0.00000182082, 6.51815E-08, -3.48553E-11, 1.11409E-14])
+air.add_model(propname='pore.molar_density',
+              model=op.models.phases.density.ideal_gas)
+
+phys = op.physics.GenericPhysics(network=pn, phase=air, geometry=geo)
+phys.add_model(propname='throat.diffusive_conductance',
+               pore_diffusivity='pore.diffusivity.pure_O2',
+               model=op.models.physics.diffusive_conductance.ordinary_diffusion)
+phys.add_model(propname='throat.hydraulic_conductance',
+               pore_diffusivity='pore.viscosity',
+               model=op.models.physics.hydraulic_conductance.classic_hagen_poiseuille)
+
+sf = op.algorithms.StokesFlow(network=pn, phase=air)
+sf.set_value_BC(pores=pn.pores('left'), values=200000)
+sf.set_value_BC(pores=pn.pores('right'), values=100000)
+sf.run()
+air.update(sf.results())
+
+air.regenerate_models()
+
+fd = op.algorithms.FickianDiffusion(network=pn, phase=air)
+fd.setup(quantity='pore.concentration.pure_O2')
+fd.set_value_BC(pores=pn.pores('left'), values=1)
+fd.set_value_BC(pores=pn.pores('right'), values=0)
+fd.run()
+air.update(fd.results())
diff --git a/openpnm/algorithms/NernstPlanck.py b/openpnm/algorithms/NernstPlanck.py
index 01a44ce30..5e82b3143 100644
--- a/openpnm/algorithms/NernstPlanck.py
+++ b/openpnm/algorithms/NernstPlanck.py
@@ -14,7 +14,7 @@ class NernstPlanck(ReactiveTransport):
                    'quantity': 'pore.concentration.'+ion,
                    'conductance': 'throat.ad_dif_mig_conductance.'+ion}
         super().__init__(**kwargs)
-        self.name = ion
+        # self.name = electrolyte  # This interfers with component name
         self.settings.update(def_set)
         self.settings.update(settings)
         if phase is not None:
diff --git a/openpnm/core/Base.py b/openpnm/core/Base.py
index b8ffdd940..b9af6f7ca 100644
--- a/openpnm/core/Base.py
+++ b/openpnm/core/Base.py
@@ -1704,16 +1704,16 @@ class Base(dict):
         return element + '.' + propname.split('.')[-1]
 
     def __str__(self):
-        horizonal_rule = '―' * 78
-        lines = [horizonal_rule]
+        horizontal_rule = '―' * 78
+        lines = [horizontal_rule]
         lines.append(self.__module__.replace('__', '') + ' : ' + self.name)
-        lines.append(horizonal_rule)
+        lines.append(horizontal_rule)
         lines.append("{0:<5s} {1:<45s} {2:<10s}".format('#',
                                                         'Properties',
                                                         'Valid Values'))
         fmt = "{0:<5d} {1:<45s} {2:>5d} / {3:<5d}"
-        lines.append(horizonal_rule)
-        props = list(set(self.keys()).difference(set(self.labels())))
+        lines.append(horizontal_rule)
+        props = self.props()
         props.sort()
         for i, item in enumerate(props):
             prop = item
@@ -1729,11 +1729,11 @@ class Base(dict):
                 defined = sp.shape(self[item])[0] \
                     - a.sum(axis=0, keepdims=(a.ndim-1) == 0)[0]
                 lines.append(fmt.format(i + 1, prop, defined, required))
-        lines.append(horizonal_rule)
+        lines.append(horizontal_rule)
         lines.append("{0:<5s} {1:<45s} {2:<10s}".format('#',
                                                         'Labels',
                                                         'Assigned Locations'))
-        lines.append(horizonal_rule)
+        lines.append(horizontal_rule)
         labels = self.labels()
         labels.sort()
         fmt = "{0:<5d} {1:<45s} {2:<10d}"
@@ -1743,7 +1743,7 @@ class Base(dict):
                 prop = prop[0:32] + '...'
             if '._' not in prop:
                 lines.append(fmt.format(i + 1, prop, sp.sum(self[item])))
-        lines.append(horizonal_rule)
+        lines.append(horizontal_rule)
         return '\n'.join(lines)
 
     def _mro(self):
diff --git a/openpnm/models/misc/__init__.py b/openpnm/models/misc/__init__.py
index 2091da78b..034ee8dc7 100644
--- a/openpnm/models/misc/__init__.py
+++ b/openpnm/models/misc/__init__.py
@@ -12,5 +12,6 @@ generating an array of random numbers.
 
 from .statistical_distributions import generic_distribution, normal, weibull, random
 from .simple_equations import linear, polynomial, generic_function
-from .basic_math import constant, product, scaled, clip, normalize
+from .basic_math import constant, product, scaled, clip, normalize, summation
+from .basic_math import fraction
 from .neighbor_lookups import from_neighbor_throats, from_neighbor_pores
diff --git a/openpnm/models/misc/basic_math.py b/openpnm/models/misc/basic_math.py
index 6f0a3df3a..ac9da674c 100644
--- a/openpnm/models/misc/basic_math.py
+++ b/openpnm/models/misc/basic_math.py
@@ -12,6 +12,49 @@ from openpnm.utils import logging
 logger = logging.getLogger(__name__)
 
 
+def fraction(target, numerator, denominator):
+    r"""
+    Calculates the ratio between two values
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    numerator : string
+        Dictionary key pointing the numerator values
+
+    denominator : string
+        Dictionary key pointing the denominator values
+    """
+    x = target[numerator]
+    y = target[denominator]
+    return x/y
+
+
+def summation(target, props=[]):
+    r"""
+    Sums the values in the given arrays
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    props : list of dictionary keys
+        The dictionary keys pointing the arrays whose elements should be summed
+
+    """
+    vals = np.zeros_like(target[props[0]])
+    for item in props:
+        vals += target[item]
+    return vals
+
+
 def normalize(target, prop, xmin=0, xmax=1):
     r"""
     Normalizes the given array between the supplied limits
diff --git a/openpnm/models/phases/__init__.py b/openpnm/models/phases/__init__.py
index c680a7c51..3894fd9e0 100644
--- a/openpnm/models/phases/__init__.py
+++ b/openpnm/models/phases/__init__.py
@@ -18,3 +18,4 @@ from . import thermal_conductivity
 from . import vapor_pressure
 from . import viscosity
 from . import misc
+from . import mixtures
diff --git a/openpnm/models/phases/density.py b/openpnm/models/phases/density.py
index d04c740ab..db45ebfe7 100644
--- a/openpnm/models/phases/density.py
+++ b/openpnm/models/phases/density.py
@@ -14,6 +14,11 @@ def standard(target, mol_weight='pore.molecular_weight',
 
     Parameters
     ----------
+    target : OpenPNM Object
+        The object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides
+        access to other necessary thermofluid properties.
+
     mol_weight : string
         The dictionary key containing the molecular weight values (kg/mol)
 
diff --git a/openpnm/models/phases/mixtures.py b/openpnm/models/phases/mixtures.py
new file mode 100644
index 000000000..7454ab08e
--- /dev/null
+++ b/openpnm/models/phases/mixtures.py
@@ -0,0 +1,200 @@
+import numpy as np
+
+
+def salinity(target, temperature='pore.temperature',
+             concentration='pore.concentration'):
+    r"""
+    Calculates the salinity in g salt per kg of solution from concentration
+
+    Parameter
+    ---------
+    target : OpenPNM Object
+        The object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides
+        access to other necessary thermofluid properties.
+    temperature : string
+        The dictionary key containing the temperature values in Kelvin.
+    concentration : string
+        The dictionary key containing the concentration values, in SI units of
+        mol/m3.
+
+    Returns
+    -------
+    salinity : ND-array
+        The salinity in g of solute per kg of solution.
+
+    Notes
+    -----
+    This model is useful for converting known concentration values (e.g.
+    calculated by a transport algorithm) into salinity values, which can then
+    be used for finding other physical properties of water which are available
+    as a function of salinity.
+
+    The salinity correlations are valid for salinity up to 160 g/kg, which
+    corresponds to a concentration of 0.05 mol/L (assuming NaCl is the only
+    solute)
+
+    """
+    C = target[concentration]
+    T = target[temperature]
+    a = 8.73220929e+00
+    b = 6.00389629e+01
+    c = -1.19083743e-01
+    d = -1.77796042e+00
+    e = 3.26987130e-04
+    f = -1.09636011e-01
+    g = -1.83933426e-07
+    S = a + b*C + c*T + d*C**2 + e*T**2 + f*C**3 + g*T**3
+    return S
+
+
+def mole_weighted_average(target, prop):
+    r"""
+
+    Parameters
+    ----------
+    target : OpenPNM Mixture Object
+        The object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides
+        access to other necessary thermofluid properties.
+    prop : string
+        The dictionary key to the property to be averaged.
+
+    Returns
+    -------
+    vals : ND-array
+        An ND-array containing the mole fraction weighted average value of the
+        specified property.
+    """
+    comps = target.components.values()
+    element = prop.split('.')[0]
+    vals = np.zeros(target._count(element))
+    for item in comps:
+        frac = target[element + '.mole_fraction.' + item.name]
+        temp = item[prop]
+        vals += temp*frac
+    return vals
+
+
+def fuller_diffusivity(target, molecular_weight='pore.molecular_weight',
+                       molar_diffusion_volume='pore.molar_diffusion_volume',
+                       temperature='pore.temperature',
+                       pressure='pore.pressure'):
+    r"""
+    Estimates the diffusion coeffient of both species in a binary gas
+    mixture using the Fuller correlation
+
+    Parameters
+    ----------
+    target : OpenPNM Mixture Object
+        The object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides
+        access to other necessary thermofluid properties.
+    molecular_weight : string
+        Dictionary key containing the molecular weight of each species.  The
+        default is 'pore.molecular_weight'
+    molar_diffusion_volume : string
+        Dictionary key containing the molar diffusion volume of each species.
+        This is used by the Fuller correlation.  The default is
+        'pore.molar_diffusion_volume'
+    temperature : string
+        Dictionary key contain the temperature of the mixture.  The default
+        is 'pore.temperature'
+    pressure : string
+        Dictionary key contain the pressure of the mixture.  The default
+        is 'pore.pressure'
+
+    Returns
+    -------
+    Dij : dict containing ND-arrys
+        The dict contains one array for each component, containing the
+        diffusion coefficient of that component at each location.
+    """
+    species_A, species_B = target.components.values()
+    T = target[temperature]
+    P = target[pressure]
+    MA = species_A[molecular_weight]
+    MB = species_B[molecular_weight]
+    vA = species_A[molar_diffusion_volume]
+    vB = species_B[molar_diffusion_volume]
+    MAB = 1e3*2*(1.0/MA + 1.0/MB)**(-1)
+    P = P*1e-5
+    value = 0.00143*T**1.75/(P*(MAB**0.5)*(vA**(1./3) + vB**(1./3))**2)*1e-4
+    return value
+
+
+def wilke_fuller_diffusivity(target, molecular_weight='pore.molecular_weight',
+                             molar_diffusion_volume='pore.molar_diffusion_volume',
+                             temperature='pore.temperature',
+                             pressure='pore.pressure'):
+    r"""
+    Estimates the diffusion coeffient of each species in a gas mixture
+
+    Uses the fuller equation to estimate binary diffusivity between pairs, then
+    uses the correction of Fairbanks and Wilke to account for the composition
+    of the gas mixture.
+
+    Parameters
+    ----------
+    target : OpenPNM Mixture Object
+        The object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides
+        access to other necessary thermofluid properties.
+    molecular_weight : string
+        Dictionary key containing the molecular weight of each species.  The
+        default is 'pore.molecular_weight'
+    molar_diffusion_volume : string
+        Dictionary key containing the molar diffusion volume of each species.
+        This is used by the Fuller correlation.  The default is
+        'pore.molar_diffusion_volume'
+    temperature : string
+        Dictionary key contain the temperature of the mixture.  The default
+        is 'pore.temperature'
+    pressure : string
+        Dictionary key contain the pressure of the mixture.  The default
+        is 'pore.pressure'
+
+    Returns
+    -------
+    Dij : dict containing ND-arrys
+        The dict contains one array for each component, containing the
+        diffusion coefficient of that component at each location.
+
+    Reference
+    ---------
+    Fairbanks DF and CR Wilke, Diffusion Coefficients in Multicomponent
+    Gas Mixtures. Industrial & Engineering Chemistry, 42(3), p471–475 (1950).
+    `DOI: 10.1021/ie50483a022 <http://doi.org/10.1021/ie50483a022>`_
+    """
+    comps = list(target.components.values())
+    values = {}
+    for i in range(len(comps)):
+        A = comps[i]
+        denom = 0.0
+        for j in range(len(comps)):
+            if i != j:
+                B = comps[j]
+                temp = MixDict(target=target, components=(A, B))
+                D = fuller_diffusivity(target=temp,
+                                       molecular_weight=molecular_weight,
+                                       molar_diffusion_volume=molar_diffusion_volume,
+                                       temperature=temperature,
+                                       pressure=pressure)
+                yB = target['pore.mole_fraction.' + B.name]
+                denom += yB/D
+        yA = target['pore.mole_fraction.' + A.name]
+        values[A.name] = (1 - yA)/denom
+    return values
+
+
+class MixDict(dict):
+    r"""
+    This utility dict is used to create a temporary mixture object containing
+    only two components of a mixture that has several.  This is necessary for
+    use of the fuller model for calculating binary diffusivities.
+    """
+    def __init__(self, target, components):
+        super().__init__(target)
+        self.components = {}
+        for item in components:
+            self.components.update({item.name: item})
diff --git a/openpnm/models/phases/molar_density.py b/openpnm/models/phases/molar_density.py
index e48b02f7f..786ba9b9d 100644
--- a/openpnm/models/phases/molar_density.py
+++ b/openpnm/models/phases/molar_density.py
@@ -21,10 +21,10 @@ def standard(target, mol_weight='pore.molecular_weight',
         controls the length of the calculated array, and also provides
         access to other necessary thermofluid properties.
 
-    pore_MW : string
+    mol_weight : string
         The dictionary key containing the molecular weight in kg/mol
 
-    pore_temperature : string
+    density : string
         The dictionary key containing the density in kg/m3
 
     Returns
diff --git a/openpnm/phases/Mixtures.py b/openpnm/phases/Mixtures.py
deleted file mode 100644
index fb0ce5b91..000000000
--- a/openpnm/phases/Mixtures.py
+++ /dev/null
@@ -1,43 +0,0 @@
-from openpnm.phases import GenericPhase
-import openpnm.models as mods
-
-
-class SalineWater(GenericPhase):
-    r"""
-    """
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
-
-        self['pore.molecular_weight'] = 0.01802
-        self['pore.critical_pressure'] = 2.2064E7
-        self['pore.critical_temperature'] = 647.1
-        self['pore.critical_volume'] = 0.003106
-        self['pore.contact_angle'] = 110.0
-        self['pore.electrical_conductivity'] = 1e-15
-
-        self['pore.permittivity'] = 78.303
-        self['pore.diffusivity.Na'] = 1.33e-09
-        self['pore.diffusivity.Cl'] = 2.03e-09
-        self['pore.valence.Na'] = 1
-        self['pore.valence.Cl'] = -1
-
-        self['throat.permittivity'] = 78.303
-        self['throat.diffusivity.Na'] = 1.33e-09
-        self['throat.diffusivity.Cl'] = 2.03e-09
-        self['throat.valence.Na'] = 1
-        self['throat.valence.Cl'] = -1
-
-        self.add_model(propname='pore.density',
-                       model=mods.phases.density.water)
-        self.add_model(propname='pore.molar_density',
-                       model=mods.phases.molar_density.standard)
-        self.add_model(propname='pore.surface_tension',
-                       model=mods.phases.surface_tension.water)
-        self.add_model(propname='pore.thermal_conductivity',
-                       model=mods.phases.thermal_conductivity.water)
-        self.add_model(propname='pore.vapor_pressure',
-                       model=mods.phases.vapor_pressure.antoine,
-                       A=8.088, B=1750.71, C=236.191)
-        self.add_model(propname='pore.viscosity',
-                       model=mods.phases.viscosity.water)
-        self.regenerate_models()
diff --git a/openpnm/phases/Multiphase.py b/openpnm/phases/Multiphase.py
index 3bd0aacbb..a26103d58 100644
--- a/openpnm/phases/Multiphase.py
+++ b/openpnm/phases/Multiphase.py
@@ -57,6 +57,9 @@ class MultiPhase(GenericPhase):
                               })
         self.settings.update(settings)
 
+        logger.warning('Multiphases are a beta feature and functionality ' +
+                       'may change in future versions')
+
         self['pore.occupancy.all'] = np.zeros(self.Np, dtype=float)
         self['throat.occupancy.all'] = np.zeros(self.Nt, dtype=float)
 
@@ -86,7 +89,8 @@ class MultiPhase(GenericPhase):
             self['throat.occupancy.all'] = np.sum(dict_, axis=0)
 
     def _get_phases(self):
-        phases = [self.project[item] for item in self.settings['phases']]
+        phases = {self.project[item].name: self.project[item]
+                  for item in self.settings['phases']}
         return phases
 
     phases = property(fget=_get_phases)
@@ -119,7 +123,7 @@ class MultiPhase(GenericPhase):
                                 element + 's')
         vals = np.zeros([self._count(element=element)], dtype=float)
         try:
-            for phase in self.phases:
+            for phase in self.phases.values():
                 vals += phase[prop]*self[element + '.occupancy.' + phase.name]
         except KeyError:
             vals = super().interleave_data(prop)
diff --git a/openpnm/phases/Water.py b/openpnm/phases/Water.py
index 8822e74a2..628071c48 100644
--- a/openpnm/phases/Water.py
+++ b/openpnm/phases/Water.py
@@ -93,6 +93,7 @@ class Water(GenericPhase):
         self['pore.contact_angle'] = 110.0
         self['pore.electrical_conductivity'] = 1e-15
         self['pore.diffusivity'] = 1e-9
+
         self.add_model(propname='pore.density',
                        model=mods.phases.density.water)
         self.add_model(propname='pore.molar_density',
diff --git a/openpnm/phases/__init__.py b/openpnm/phases/__init__.py
index 01cd79e61..4aade75bb 100644
--- a/openpnm/phases/__init__.py
+++ b/openpnm/phases/__init__.py
@@ -81,4 +81,3 @@ from .Air import Air
 from .Water import Water
 from .Mercury import Mercury
 from .Multiphase import MultiPhase
-from .Mixtures import SalineWater
diff --git a/openpnm/phases/mixtures/DryAir.py b/openpnm/phases/mixtures/DryAir.py
new file mode 100644
index 000000000..ea3a1e95f
--- /dev/null
+++ b/openpnm/phases/mixtures/DryAir.py
@@ -0,0 +1,25 @@
+from openpnm.phases.mixtures import GenericMixture, species
+from openpnm.models.phases.mixtures import mole_weighted_average
+import openpnm.models as mods
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+class DryAir(GenericMixture):
+    r"""
+    """
+    def __init__(self, network, **kwargs):
+        super().__init__(network=network, components=[], **kwargs)
+
+        N2 = species.gases.N2(network=network, name='N2_'+self.name)
+        O2 = species.gases.O2(network=network, name='O2_'+self.name)
+        self.settings['components'] = [O2.name, N2.name]
+        self.set_mole_fraction(component=N2, values=0.791)
+        self.set_mole_fraction(component=O2, values=0.209)
+        self.add_model(propname='pore.diffusivity.N2',
+                       model=mods.phases.mixtures.fuller_diffusivity)
+        self.add_model(propname='pore.diffusivity.O2',
+                       model=mods.phases.mixtures.fuller_diffusivity)
+        self.add_model(propname='pore.molar_mass',
+                       model=mole_weighted_average,
+                       prop='pore.molecular_weight')
diff --git a/openpnm/phases/mixtures/GenericMixture.py b/openpnm/phases/mixtures/GenericMixture.py
new file mode 100644
index 000000000..b223ccf42
--- /dev/null
+++ b/openpnm/phases/mixtures/GenericMixture.py
@@ -0,0 +1,315 @@
+# from collections import ChainMap  # Might use eventually
+import numpy as np
+from openpnm.phases import GenericPhase as GenericPhase
+from openpnm.utils import logging, HealthDict, PrintableList
+logger = logging.getLogger(__name__)
+
+
+class GenericMixture(GenericPhase):
+    r"""
+    Creates Phase object that represents a multicomponent mixture system
+    consisting of a given list of OpenPNM Phase objects as components.
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    components : list of OpenPNM Phase objects
+        A list of all components that constitute this mixture
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project
+        is used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.
+        If no name is given, one is generated.
+
+    """
+    def __init__(self, components=[], settings={}, **kwargs):
+        self.settings.update({'components': [],
+                              })
+        super().__init__(settings={'prefix': 'mix'}, **kwargs)
+        self.settings.update(settings)
+
+        # Add any supplied phases to the phases list
+        for comp in components:
+            self.settings['components'].append(comp.name)
+            self['pore.mole_fraction.'+comp.name] = 0.0
+
+        self['pore.mole_fraction.all'] = np.zeros(self.Np, dtype=float)
+        logger.warning('Mixtures are a beta feature and functionality may ' +
+                       'change in future versions')
+
+    def __getitem__(self, key):
+        try:
+            vals = super().__getitem__(key)
+        except KeyError:
+            try:
+                # If key ends in component name, fetch it
+                if key.split('.')[-1] in self.settings['components']:
+                    comp = self.project[key.split('.')[-1]]
+                    vals = comp[key.rsplit('.', maxsplit=1)[0]]
+                    return vals
+                else:
+                    raise KeyError
+            except KeyError:
+                vals = self.interleave_data(key)
+        return vals
+
+    def __setitem__(self, key, value):
+        # Prevent writing 'element.property.component' on mixture
+        invalid_keys = set(self.props(deep=True)).difference(set(self.props()))
+        if key in invalid_keys:
+            raise Exception(key + ' already assigned to a component object')
+        super().__setitem__(key, value)
+
+    def props(self, deep=False, **kwargs):
+        temp = PrintableList()
+        if deep:
+            for item in self.components.values():
+                temp.extend([prop+'.'+item.name for prop in item.props(**kwargs)])
+        temp.extend(super().props(**kwargs))
+        temp.sort()
+        return temp
+
+    def __str__(self):
+        horizontal_rule = '―' * 78
+        lines = super().__str__()
+        lines = '\n'.join((lines, 'Component Phases', horizontal_rule))
+        for item in self.components.values():
+            lines = '\n'.join((lines, item.__module__.replace('__', '') +
+                               ' : ' + item.name))
+        lines = '\n'.join((lines, horizontal_rule))
+        return lines
+
+    def _update_total_molfrac(self):
+        # Update mole_fraction.all
+        self['pore.mole_fraction.all'] = 0.0
+        dict_ = list(self['pore.mole_fraction'].values())
+        if len(dict_) > 1:
+            self['pore.mole_fraction.all'] = np.sum(dict_, axis=0)
+        self['throat.mole_fraction.all'] = 0.0
+        dict_ = list(self['throat.mole_fraction'].values())
+        if len(dict_) > 1:
+            self['throat.mole_fraction.all'] = np.sum(dict_, axis=0)
+
+    def update_concentrations(self, mole_fraction='pore.mole_fraction'):
+        r"""
+        Re-calculates the concentration of each species in the mixture based
+        on the current mole fractions.
+
+        This method looks up the mole fractions *and* the density of the
+        mixture, then finds the respective concentrations in $mol/m^{3}$.
+
+        Parameters
+        ----------
+
+        """
+        density = self['pore.molar_density']
+        for item in self.components.values():
+            mf = self['pore.mole_fraction.'+item.name]
+            self['pore.concentration.'+item.name] = density*mf
+
+    def update_mole_fractions(self, concentration='pore.concentration',
+                              molar_density='pore.molar_density'):
+        r"""
+        Re-calculates mole fraction of each species in mixture based on the
+        current concentrations.
+
+        This method looks up the concentration of each species (using the
+        optionally specified concentration dictionary key), and calculates
+        the mole fraction.  Optionally, it can use a molar density for the
+        mixture and N-1 concentrations to determine the Nth concentration and
+        all species mole fractions.
+
+        Parameters
+        ----------
+        concentration : string, optional
+            The dictionary key pointing to the desired concentration values.
+            The default is 'pore.concentration'.  Given this value, lookups
+            are performed for each species in the mixture.
+        molar_density : string, optional
+            The dictionary key pointing to the molar density of the mixture.
+            If not given (default), then 'pore.molar_density' is used. If
+            there are N-1 concentrations specified, then ``molar_density`` is
+            automatically used to find the Nth concentration.
+
+        Notes
+        -----
+        The method does not return any values.  Instead it updates the mole
+        fraction arrays of each species directly.
+        """
+        concentrations = [concentration + '.' + comp for comp
+                          in self.settings['components']
+                          if concentration + '.' + comp in self.keys()]
+        if len(concentrations) == len(self.components):
+            # Find total number of moles per unit volume
+            density = 0.0
+            for conc in concentrations:
+                density += self[conc]
+            # Normalize moles per unit volume for each species by the total
+            for conc in concentrations:
+                element, quantity, component = conc.split('.')
+                self[element+'.mole_fraction.'+component] = self[conc]/density
+        elif len(concentrations) == (len(self.components) - 1):
+            # Find mole fraction of N-1 species
+            mol_frac = 0.0
+            density = self[molar_density]
+            for conc in concentrations:
+                element, quantity, component = conc.split('.')
+                self[element+'.mole_fraction.'+component] = self[conc]/density
+                mol_frac += self[element+'.mole_fraction.'+component]
+            # Find mole fraction of Nth species using molar_density
+            given_comps = [conc.split('.')[2] for conc in concentrations]
+            all_comps = self.settings['components']
+            component = list(set(all_comps).difference(set(given_comps)))[0]
+            self[element+'.mole_fraction.'+component] = 1 - mol_frac
+            # [self[element+'.concentration.'+component] = (1 - mol_frac)*density
+        else:
+            raise Exception('Insufficient concentration values found ' +
+                            'for component species, must specify ' +
+                            str(abs(n_spec + 1)) + ' additional values')
+
+    def set_concentration(self, component, values=[]):
+        r"""
+        Specify mole fraction of each component in each pore
+
+        Parameters
+        ----------
+        components : OpenPNM Phase object or name string
+            The phase whose mole fraction is being specified
+
+        values : array_like
+            The concentration of the given ``component `` in each pore.  This
+            array must be *Np*-long, with one value for each pore in the
+            network.  If a scalar is received it is applied to all pores.
+
+        See Also
+        --------
+        set_mole_fraction
+
+        """
+        if type(component) == str:
+            component = self.components[component]
+        Pvals = np.array(values, ndmin=1)
+        if component not in self.project:
+            raise Exception(f"{component.name} doesn't belong to this project")
+        else:
+            if component.name not in self.settings['components']:
+                self.settings['components'].append(component.name)
+        if np.any(Pvals < 0.0):
+            logger.warning('Received values contain negative concentrations')
+        if Pvals.size:
+            self['pore.concentration.' + component.name] = Pvals
+
+    def set_mole_fraction(self, component, values=[]):
+        r"""
+        Specify mole fraction of each component in each pore
+
+        Parameters
+        ----------
+        components : OpenPNM Phase object or name string
+            The phase whose mole fraction is being specified
+
+        values : array_like
+            The mole fraction of the given ``component `` in each pore.  This
+            array must be *Np*-long, with one value between 0 and 1 for each
+            pore in the network.  If a scalar is received it is applied to
+            all pores.
+
+        See Also
+        --------
+        set_concentration
+
+        """
+        if type(component) == str:
+            component = self.components[component]
+        Pvals = np.array(values, ndmin=1)
+        if component not in self.project:
+            raise Exception(f"{component.name} doesn't belong to this project")
+        else:
+            if component.name not in self.settings['components']:
+                self.settings['components'].append(component.name)
+        if np.any(Pvals > 1.0) or np.any(Pvals < 0.0):
+            logger.warning('Received values contain mole fractions outside ' +
+                           'the range of 0 to 1')
+        if Pvals.size:
+            self['pore.mole_fraction.' + component.name] = Pvals
+        self._update_total_molfrac()
+
+    def _get_comps(self):
+        comps = {item: self.project[item] for item in self.settings['components']}
+        return comps
+
+    def _set_comps(self, components):
+        if not isinstance(components, list):
+            components = [components]
+        self.settings['components'] = [val.name for val in components]
+
+    components = property(fget=_get_comps, fset=_set_comps)
+
+    def interleave_data(self, prop):
+        r"""
+        Gathers property values from component phases to build a single array
+
+        If the requested ``prop`` is not on this Mixture, then a search is
+        conducted on all associated components objects, and values from each
+        are assembled into a single array.
+
+        Parameters
+        ----------
+        prop : string
+            The property to be retrieved
+
+        Returns
+        -------
+        array : ND-array
+            An array containing the specified property retrieved from each
+            component phase and assembled based on the specified mixing rule
+
+        """
+        element = prop.split('.')[0]
+        if element == 'pore':
+            if np.any(self[element + '.mole_fraction.all'] != 1.0):
+                self._update_total_molfrac()
+                if np.any(self[element + '.mole_fraction.all'] != 1.0):
+                    raise Exception('Mole fraction does not add to unity in all ' +
+                                    element + 's')
+        vals = np.zeros([self._count(element=element)], dtype=float)
+        try:
+            for comp in self.components.values():
+                vals += comp[prop]*self[element+'.mole_fraction.'+comp.name]
+        except KeyError:
+            vals = super().interleave_data(prop)
+        return vals
+
+    def check_mixture_health(self):
+        r"""
+        Checks the "health" of the mixture
+
+        Calculates the mole fraction of all species in each pore and returns
+        an list of where values are too low or too high
+
+        Returns
+        -------
+        health : dict
+            A HealtDict object containing lists of locations where the mole
+            fractions are not unity.  One value indiates locations that are
+            too high, and another where they are too low.
+
+        """
+        h = HealthDict()
+        h['mole_fraction_too_low'] = []
+        h['mole_fraction_too_high'] = []
+        self._update_total_molfrac()
+        lo = np.where(self['pore.mole_fraction.all'] < 1.0)[0]
+        hi = np.where(self['pore.mole_fraction.all'] > 1.0)[0]
+        if len(lo) > 0:
+            h['mole_fraction_too_low'] = lo
+        if len(hi) > 0:
+            h['mole_fraction_too_high'] = hi
+        return h
diff --git a/openpnm/phases/mixtures/GenericSpecies.py b/openpnm/phases/mixtures/GenericSpecies.py
new file mode 100644
index 000000000..fe99c14b9
--- /dev/null
+++ b/openpnm/phases/mixtures/GenericSpecies.py
@@ -0,0 +1,33 @@
+from openpnm.phases import GenericPhase as GenericPhase
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+class GenericSpecies(GenericPhase):
+    r"""
+    Creates Phase object that represents a single species in a multicomponent
+    mixture system.
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project
+        is used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.
+        If no name is given, one is generated.
+
+    """
+    @property
+    def mixture(self):
+        for item in self.project:
+            if hasattr(item, 'components'):
+                for comp in item.components.values():
+                    if self is comp:
+                        return item
diff --git a/openpnm/phases/mixtures/HumidAir.py b/openpnm/phases/mixtures/HumidAir.py
new file mode 100644
index 000000000..f6de1bf6d
--- /dev/null
+++ b/openpnm/phases/mixtures/HumidAir.py
@@ -0,0 +1,30 @@
+from openpnm.phases.mixtures import GenericMixture, species
+import openpnm.models as mods
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+class HumidAir(GenericMixture):
+    r"""
+    """
+    def __init__(self, network, **kwargs):
+        super().__init__(network=network, components=[], **kwargs)
+
+        N2 = species.gases.N2(network=network, name='N2_'+self.name)
+        O2 = species.gases.O2(network=network, name='O2_'+self.name)
+        H2O = species.liquids.H2O(network=network, name='H2O_'+self.name)
+        self.settings['components'] = [O2.name, N2.name, H2O.name]
+        self.set_mole_fraction(component=N2, values=0.791)
+        self.set_mole_fraction(component=O2, values=0.209)
+        self.set_mole_fraction(component=H2O, values=0.000)
+        self.add_model(propname='pore.molar_mass',
+                       model=mods.phases.mixtures.mole_weighted_average,
+                       prop='pore.molecular_weight')
+        self.add_model(propname='pore.vapor_pressure',
+                       model=mods.phases.vapor_pressure.water)
+#        self.add_model(propname='pore.mole_fraction.'+H2O.name,
+#                       model=mods.misc.fraction,
+#                       numerator='pore.vapor_pressure',
+#                       denominator='pore.pressure')
+        self.add_model(propname='pore.diffusivity',
+                       model=mods.phases.mixtures.wilke_fuller_diffusivity)
diff --git a/openpnm/phases/mixtures/SalineWater.py b/openpnm/phases/mixtures/SalineWater.py
new file mode 100644
index 000000000..5ea34efe0
--- /dev/null
+++ b/openpnm/phases/mixtures/SalineWater.py
@@ -0,0 +1,35 @@
+from openpnm.phases.mixtures import GenericMixture, species
+from openpnm import models as mods
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+class SalineWater(GenericMixture):
+    r"""
+    """
+    def __init__(self, network, components=None, **kwargs):
+        if components is not None:
+            logger.warn('Ignoring received components')
+        super().__init__(network=network, components=[], **kwargs)
+
+        Cl = species.ions.Cl(network=network, name='Cl_'+self.name)
+        Na = species.ions.Na(network=network, name='Na_'+self.name)
+        W = species.liquids.H2O(network=network, name='H2O_'+self.name)
+        self.settings['components'] = [Cl.name, Na.name, W.name]
+        self.set_concentration(component=W, values=998/0.018)
+        self.set_concentration(component=Na, values=0.0)
+        self.set_concentration(component=Cl, values=0.0)
+        self.add_model(propname='pore.salt_concentration',
+                       model=mods.misc.summation,
+                       props=['pore.concentration.Na_'+self.name,
+                              'pore.concentration.Cl_'+self.name])
+        self.add_model(propname='pore.salinity',
+                       model=mods.phases.mixtures.salinity,
+                       concentration='pore.concentration.Na_'+self.name)
+        self.add_model(propname='pore.mass_density',
+                       model=mods.phases.density.water,
+                       salinity='pore.salinity')
+        self.add_model(propname='pore.viscosity',
+                       model=mods.phases.viscosity.water,
+                       salinity='pore.salinity')
+        self.update_mole_fractions()
diff --git a/openpnm/phases/mixtures/__init__.py b/openpnm/phases/mixtures/__init__.py
new file mode 100644
index 000000000..c45150023
--- /dev/null
+++ b/openpnm/phases/mixtures/__init__.py
@@ -0,0 +1,6 @@
+from .GenericMixture import GenericMixture
+from .GenericSpecies import GenericSpecies
+from . import species
+from .SalineWater import SalineWater
+from .DryAir import DryAir
+from .HumidAir import HumidAir
diff --git a/openpnm/phases/mixtures/species/__init__.py b/openpnm/phases/mixtures/species/__init__.py
new file mode 100644
index 000000000..75525186f
--- /dev/null
+++ b/openpnm/phases/mixtures/species/__init__.py
@@ -0,0 +1,3 @@
+from . import ions
+from . import gases
+from . import liquids
diff --git a/openpnm/phases/mixtures/species/gases/CO2.py b/openpnm/phases/mixtures/species/gases/CO2.py
new file mode 100644
index 000000000..8894c4ebc
--- /dev/null
+++ b/openpnm/phases/mixtures/species/gases/CO2.py
@@ -0,0 +1,36 @@
+from openpnm.phases import GenericPhase
+import openpnm.models as mods
+
+
+class CO2(GenericPhase):
+    r"""
+    Creates Phase object with preset models and values for CO2 gas
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> CO2 = mixtures.species.gases.CO2(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self['pore.molecular_weight'] = 0.04401  # kg/mol
+        self['pore.molar_diffusion_volume'] = 17.9  # Wrong
diff --git a/openpnm/phases/mixtures/species/gases/H2.py b/openpnm/phases/mixtures/species/gases/H2.py
new file mode 100644
index 000000000..8d6e6a4c1
--- /dev/null
+++ b/openpnm/phases/mixtures/species/gases/H2.py
@@ -0,0 +1,36 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class H2(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for Na ions
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> H2 = mixtures.species.gases.H2(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self['pore.molecular_weight'] = 0.002016  # kg/mol
+        self['pore.molar_diffusion_volume'] = 17.9  # Wrong
diff --git a/openpnm/phases/mixtures/species/gases/N2.py b/openpnm/phases/mixtures/species/gases/N2.py
new file mode 100644
index 000000000..dfa6be18a
--- /dev/null
+++ b/openpnm/phases/mixtures/species/gases/N2.py
@@ -0,0 +1,38 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class N2(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for N2 gas
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> N2 = mixtures.species.gases.N2(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.pop('pore.temperature', None)
+        self.pop('pore.pressure', None)
+
+        self['pore.molecular_weight'] = 0.0280  # kg/mol
+        self['pore.molar_diffusion_volume'] = 17.9  # ??
diff --git a/openpnm/phases/mixtures/species/gases/O2.py b/openpnm/phases/mixtures/species/gases/O2.py
new file mode 100644
index 000000000..ca5a32d39
--- /dev/null
+++ b/openpnm/phases/mixtures/species/gases/O2.py
@@ -0,0 +1,38 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class O2(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for O2 gas
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> O2 = mixtures.species.gases.O2(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.pop('pore.temperature', None)
+        self.pop('pore.pressure', None)
+
+        self['pore.molecular_weight'] = 0.032  # kg/mol
+        self['pore.molar_diffusion_volume'] = 16.6  # ??
diff --git a/openpnm/phases/mixtures/species/gases/__init__.py b/openpnm/phases/mixtures/species/gases/__init__.py
new file mode 100644
index 000000000..11bf827d0
--- /dev/null
+++ b/openpnm/phases/mixtures/species/gases/__init__.py
@@ -0,0 +1,4 @@
+from .N2 import N2
+from .O2 import O2
+from .CO2 import CO2
+from .H2 import H2
diff --git a/openpnm/phases/mixtures/species/ions/Cl.py b/openpnm/phases/mixtures/species/ions/Cl.py
new file mode 100644
index 000000000..605252ccf
--- /dev/null
+++ b/openpnm/phases/mixtures/species/ions/Cl.py
@@ -0,0 +1,39 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class Cl(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for Cl ions
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> Cl = mixtures.species.ions.Cl(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self['pore.molecular_weight'] = 0.03545  # kg/mol
+        self['pore.diffusivity'] = 2.03e-09  # m2/s
+        self['throat.diffusivity'] = 2.03e-09  # m2/s
+        self['pore.valence'] = -1
+        self['throat.valence'] = -1
diff --git a/openpnm/phases/mixtures/species/ions/Na.py b/openpnm/phases/mixtures/species/ions/Na.py
new file mode 100644
index 000000000..8140b6692
--- /dev/null
+++ b/openpnm/phases/mixtures/species/ions/Na.py
@@ -0,0 +1,39 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class Na(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for Na ions
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> Na = mixtures.species.ions.Na(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self['pore.molecular_weight'] = 0.022990  # kg/mol
+        self['pore.diffusivity'] = 1.33e-09  # m2/s
+        self['throat.diffusivity'] = 1.33e-09  # m2/s
+        self['pore.valence'] = 1
+        self['throat.valence'] = 1
diff --git a/openpnm/phases/mixtures/species/ions/Proton.py b/openpnm/phases/mixtures/species/ions/Proton.py
new file mode 100644
index 000000000..5ef0426a8
--- /dev/null
+++ b/openpnm/phases/mixtures/species/ions/Proton.py
@@ -0,0 +1,36 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class Proton(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for H^+ ions
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> H = mixtures.species.ions.Proton(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self['pore.molecular_weight'] = 0.0291
+        self['pore.diffusivity'] = 0.1
diff --git a/openpnm/phases/mixtures/species/ions/__init__.py b/openpnm/phases/mixtures/species/ions/__init__.py
new file mode 100644
index 000000000..44c70a7f5
--- /dev/null
+++ b/openpnm/phases/mixtures/species/ions/__init__.py
@@ -0,0 +1,3 @@
+from .Cl import Cl
+from .Na import Na
+from .Proton import Proton
diff --git a/openpnm/phases/mixtures/species/liquids/H2O.py b/openpnm/phases/mixtures/species/liquids/H2O.py
new file mode 100644
index 000000000..458940fee
--- /dev/null
+++ b/openpnm/phases/mixtures/species/liquids/H2O.py
@@ -0,0 +1,36 @@
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class H2O(mixtures.GenericSpecies):
+    r"""
+    Creates Phase object with preset models and values for H2O ions
+
+    Parameters
+    ----------
+    network : OpenPNM Network object
+        The network to which this phase object will be attached.
+
+    project : OpenPNM Project object, optional
+        The Project with which this phase should be associted.  If a
+        ``network`` is given then this is ignored and the Network's project is
+        used.  If a ``network`` is not given then this is mandatory.
+
+    name : string, optional
+        The name of the phase.  This is useful to keep track of the objects
+        throughout the simulation.  The name must be unique to the project.  If
+        no name is given, one is generated.
+
+    Examples
+    --------
+    >>> import openpnm as op
+    >>> import openpnm.phases.mixtures as mixtures
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> H2O = mixtures.species.liquids.H2O(network=pn)
+
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self['pore.molecular_weight'] = 0.0291
+        self['pore.molar_diffusion_volume'] = 17.9  # Wrong
diff --git a/openpnm/phases/mixtures/species/liquids/__init__.py b/openpnm/phases/mixtures/species/liquids/__init__.py
new file mode 100644
index 000000000..19489cbab
--- /dev/null
+++ b/openpnm/phases/mixtures/species/liquids/__init__.py
@@ -0,0 +1,1 @@
+from .H2O import H2O
diff --git a/openpnm/utils/misc.py b/openpnm/utils/misc.py
index f62da57e1..df79f7edb 100644
--- a/openpnm/utils/misc.py
+++ b/openpnm/utils/misc.py
@@ -183,8 +183,12 @@ class HealthDict(PrintableDict):
     def _get_health(self):
         health = True
         for item in list(self.keys()):
-            if self[item] != []:
-                health = False
+            try:
+                if len(self[item]) > 0:
+                    health = False
+            except TypeError:
+                if self[item]:
+                    health = False
         return health
 
     health = property(fget=_get_health)
diff --git a/setup.py b/setup.py
index 79c012ead..e5fb2f10b 100644
--- a/setup.py
+++ b/setup.py
@@ -32,6 +32,10 @@ setup(
         'openpnm.network',
         'openpnm.geometry',
         'openpnm.phases',
+        'openpnm.phases.mixtures',
+        'openpnm.phases.mixtures.species.ions',
+        'openpnm.phases.mixtures.species.liquids',
+        'openpnm.phases.mixtures.species.gases',
         'openpnm.physics',
         'openpnm.utils',
         'openpnm.io',


Test patch:
diff --git a/tests/unit/models/phase/MixturesTest.py b/tests/unit/models/phase/MixturesTest.py
new file mode 100644
index 000000000..619ff84ac
--- /dev/null
+++ b/tests/unit/models/phase/MixturesTest.py
@@ -0,0 +1,70 @@
+import openpnm as op
+import numpy as np
+from openpnm.phases import mixtures
+import openpnm.models as mods
+
+
+class MixturesTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.N2 = mixtures.species.gases.N2(network=self.net)
+        self.O2 = mixtures.species.gases.O2(network=self.net)
+        self.H2O = mixtures.species.liquids.H2O(network=self.net)
+        self.ha = mixtures.GenericMixture(network=self.net,
+                                          components=[self.N2,
+                                                      self.O2,
+                                                      self.H2O])
+
+    def test_wilke_fuller_diffusivity(self):
+        self.ha.set_concentration(component=self.N2, values=0.7)
+        self.ha.set_concentration(component=self.O2, values=0.2)
+        self.ha.set_concentration(component=self.H2O, values=0.1)
+        self.ha.update_mole_fractions()
+        self.ha.add_model(propname='pore.diffusivity',
+                          model=mods.phases.mixtures.wilke_fuller_diffusivity)
+        assert 'pore.diffusivity.' + self.O2.name in self.ha.keys()
+        assert 'pore.diffusivity.' + self.N2.name in self.ha.keys()
+        assert 'pore.diffusivity.' + self.H2O.name in self.ha.keys()
+        Da = self.ha['pore.diffusivity.' + self.O2.name][0]
+        Db = self.ha['pore.diffusivity.' + self.N2.name][0]
+        Dc = self.ha['pore.diffusivity.' + self.H2O.name][0]
+        assert Da != Db
+        assert Db != Dc
+
+        self.ha.set_concentration(component=self.N2, values=0.8)
+        self.ha.set_concentration(component=self.O2, values=0.2)
+        self.ha.set_concentration(component=self.H2O, values=0.0)
+        self.ha.update_mole_fractions()
+        self.ha.regenerate_models()
+        Da = self.ha['pore.diffusivity.' + self.O2.name][0]
+        Db = self.ha['pore.diffusivity.' + self.N2.name][0]
+        Dc = self.ha['pore.diffusivity.' + self.H2O.name][0]
+        assert np.allclose(Da, Db)
+        assert Db != Dc
+
+    def test_mole_weighted_average(self):
+        self.ha.set_concentration(component=self.N2, values=0.5)
+        self.ha.set_concentration(component=self.O2, values=0.5)
+        self.ha.set_concentration(component=self.H2O, values=0.0)
+        self.N2['pore.test'] = 3
+        self.O2['pore.test'] = 3
+        self.ha.add_model(propname='pore.blah',
+                          model=mods.phases.mixtures.mole_weighted_average,
+                          prop='pore.test')
+        assert 'pore.blah' not in self.ha.keys()
+        self.H2O['pore.test'] = 3
+        self.ha.add_model(propname='pore.blah',
+                          model=mods.phases.mixtures.mole_weighted_average,
+                          prop='pore.test')
+        assert 'pore.blah' in self.ha.keys()
+
+
+if __name__ == '__main__':
+
+    t = MixturesTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/phases/mixtures/DryAirTest.py b/tests/unit/phases/mixtures/DryAirTest.py
new file mode 100644
index 000000000..8cc202f86
--- /dev/null
+++ b/tests/unit/phases/mixtures/DryAirTest.py
@@ -0,0 +1,24 @@
+import openpnm as op
+from openpnm.phases import mixtures
+
+
+class DryAirTest:
+    def setup_class(self):
+        ws = op.Workspace()
+        ws.clear()
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+
+    def test_init(self):
+        da = mixtures.DryAir(network=self.net)
+        assert isinstance(da, mixtures.GenericMixture)
+
+
+if __name__ == '__main__':
+
+    t = DryAirTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/phases/mixtures/GenericMixtureTest.py b/tests/unit/phases/mixtures/GenericMixtureTest.py
new file mode 100644
index 000000000..ff6464fde
--- /dev/null
+++ b/tests/unit/phases/mixtures/GenericMixtureTest.py
@@ -0,0 +1,106 @@
+import openpnm as op
+from openpnm.phases import mixtures
+import scipy as sp
+import pytest
+
+
+class MixtureTest:
+    def setup_class(self):
+        ws = op.Workspace()
+        ws.clear()
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+        self.N2 = mixtures.species.gases.N2(network=self.net, name='pure_N2')
+        self.O2 = mixtures.species.gases.O2(network=self.net, name='pure_O2')
+        self.CO2 = mixtures.species.gases.CO2(network=self.net, name='pure_CO2')
+        self.H2 = mixtures.species.gases.H2(network=self.net, name='pure_H2')
+        self.air = mixtures.GenericMixture(network=self.net,
+                                           components=[self.N2, self.O2,
+                                                       self.H2, self.CO2],
+                                           name='air_mixture')
+
+    def test_set_mole_fraction(self):
+        self.air.set_mole_fraction(self.N2, 0.790)
+        self.air.set_mole_fraction(self.O2, 0.209)
+        self.air.set_mole_fraction(self.CO2, 0.001)
+        assert sp.all(self.air['pore.mole_fraction.all'] == 1.0)
+
+    def test_props(self):
+        a = self.air.props(deep=False)
+        b = self.air.props(deep=True)
+        assert len(b) > len(a)
+
+    def test_update_mole_fraction_with_molar_density(self):
+        self.air.pop('pore.concentration.'+self.N2.name, None)
+        self.air['pore.concentration.'+self.O2.name] = 0.5
+        self.air['pore.concentration.'+self.CO2.name] = 0.0
+        self.air['pore.concentration.'+self.H2.name] = 0.0
+        self.air['pore.molar_density'] = 2.0
+        self.air.update_mole_fractions(molar_density='pore.molar_density')
+        assert sp.all(self.air['pore.mole_fraction.all'] == 1.0)
+
+    def test_update_mole_fraction_with_all_concentrations(self):
+        self.air['pore.concentration.'+self.O2.name] = 1.5
+        self.air['pore.concentration.'+self.N2.name] = 0.5
+        self.air['pore.concentration.'+self.CO2.name] = 0.0
+        self.air['pore.concentration.'+self.H2.name] = 0.0
+        self.air.update_mole_fractions()
+        assert sp.all(self.air['pore.mole_fraction.all'] == 1.0)
+
+    def test_interleave_data(self):
+        r"""
+        """
+        self.air['pore.concentration.'+self.O2.name] = 1.0
+        self.air['pore.concentration.'+self.N2.name] = 0.0
+        self.air['pore.concentration.'+self.CO2.name] = 0.0
+        self.air['pore.concentration.'+self.H2.name] = 0.0
+        self.air.update_mole_fractions()
+        MW = self.air['pore.molecular_weight'][0]
+        assert MW == 0.032
+        self.air['pore.concentration.'+self.O2.name] = 1.0
+        self.air['pore.concentration.'+self.N2.name] = 0.5
+        self.air['pore.concentration.'+self.CO2.name] = 0.3
+        self.air['pore.concentration.'+self.H2.name] = 0.2
+        self.air.update_mole_fractions()
+        MW = self.air['pore.molecular_weight'][0]
+        assert MW == 0.0298031
+
+    def test_check_health(self):
+        self.air.set_mole_fraction(self.N2, 0.790)
+        self.air.set_mole_fraction(self.O2, 0.209)
+        self.air.set_mole_fraction(self.CO2, 0.001)
+        self.air.set_mole_fraction(self.H2, 0.000)
+        h = self.air.check_mixture_health()
+        assert h.health is True
+        self.air.set_mole_fraction(self.CO2, 0.002)
+        h = self.air.check_mixture_health()
+        assert h.health is False
+        self.air.set_mole_fraction(self.CO2, 0.000)
+        h = self.air.check_mixture_health()
+        assert h.health is False
+        assert len(h['mole_fraction_too_low']) == self.air.Np
+        self.air.set_mole_fraction(self.CO2, 0.001)
+        self.air['pore.mole_fraction.'+self.CO2.name][0] = 0.0
+        h = self.air.check_mixture_health()
+        assert h.health is False
+        assert len(h['mole_fraction_too_low']) == 1
+
+    def test_getitem(self):
+        d = self.air['pore.mole_fraction']
+        set_a = set(['pore.mole_fraction.pure_N2',
+                     'pore.mole_fraction.pure_O2',
+                     'pore.mole_fraction.pure_H2',
+                     'pore.mole_fraction.pure_CO2',
+                     'pore.mole_fraction.all'])
+        assert set_a.difference(set(d.keys())) == set()
+
+
+if __name__ == '__main__':
+
+    t = MixtureTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/phases/mixtures/GenericSpeciesTest.py b/tests/unit/phases/mixtures/GenericSpeciesTest.py
new file mode 100644
index 000000000..31a7104a4
--- /dev/null
+++ b/tests/unit/phases/mixtures/GenericSpeciesTest.py
@@ -0,0 +1,32 @@
+import openpnm as op
+from openpnm.phases import mixtures
+
+
+class SpeciesTest:
+    def setup_class(self):
+        ws = op.Workspace()
+        ws.clear()
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+        self.N2 = mixtures.species.gases.N2(network=self.net, name='pure_N2')
+        self.O2 = mixtures.species.gases.O2(network=self.net, name='pure_O2')
+        self.CO2 = mixtures.species.gases.CO2(network=self.net, name='pure_CO2')
+        self.H2 = mixtures.species.gases.H2(network=self.net, name='pure_H2')
+        self.air = mixtures.GenericMixture(network=self.net,
+                                           components=[self.N2, self.O2,
+                                                       self.H2, self.CO2],
+                                           name='air_mixture')
+
+    def test_lookup_mixture(self):
+        assert self.N2.mixture is self.air
+
+
+if __name__ == '__main__':
+
+    t = SpeciesTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/phases/mixtures/HumidAirTest.py b/tests/unit/phases/mixtures/HumidAirTest.py
new file mode 100644
index 000000000..aca2dc478
--- /dev/null
+++ b/tests/unit/phases/mixtures/HumidAirTest.py
@@ -0,0 +1,24 @@
+import openpnm as op
+from openpnm.phases import mixtures
+
+
+class HumidAirTest:
+    def setup_class(self):
+        ws = op.Workspace()
+        ws.clear()
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+
+    def test_init(self):
+        ha = mixtures.GenericMixture(network=self.net)
+        assert isinstance(ha, mixtures.GenericMixture)
+
+
+if __name__ == '__main__':
+
+    t = HumidAirTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/phases/mixtures/SalineWaterTest.py b/tests/unit/phases/mixtures/SalineWaterTest.py
new file mode 100644
index 000000000..3ecec2cd3
--- /dev/null
+++ b/tests/unit/phases/mixtures/SalineWaterTest.py
@@ -0,0 +1,24 @@
+import openpnm as op
+from openpnm.phases import mixtures
+
+
+class SalineWaterTest:
+    def setup_class(self):
+        ws = op.Workspace()
+        ws.clear()
+        self.net = op.network.Cubic(shape=[10, 10, 10])
+
+    def test_init(self):
+        self.sw = mixtures.SalineWater(network=self.net)
+        assert isinstance(self.sw, mixtures.GenericMixture)
+
+
+if __name__ == '__main__':
+
+    t = SalineWaterTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()