Repository: pybamm-team/PyBaMM
Pull request #525 (2019-07-09 15:08:47)

Issue 490 side reactions
# Description

Add side reactions for leading-order and full lead-acid models.
Fixes some of #490 

## Type of change

- [x] New feature (non-breaking change which adds functionality)

# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks: 

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works
- [x] Any dependent changes have been merged and published in downstream modules



Patch:
diff --git a/docs/source/models/submodels/index.rst b/docs/source/models/submodels/index.rst
index 5188eb479..ff2cc824d 100644
--- a/docs/source/models/submodels/index.rst
+++ b/docs/source/models/submodels/index.rst
@@ -5,11 +5,12 @@ Submodels
     :maxdepth: 1
 
     base_submodel
+    current_collector/index
+    convection/index
     electrode/index
     electrolyte/index
     interface/index
+    oxygen_diffusion/index
     particle/index
     porosity/index
     thermal/index
-    current_collector/index
-    convection/index
diff --git a/docs/source/models/submodels/oxygen_diffusion/base_oxygen_diffusion.rst b/docs/source/models/submodels/oxygen_diffusion/base_oxygen_diffusion.rst
new file mode 100644
index 000000000..fef573564
--- /dev/null
+++ b/docs/source/models/submodels/oxygen_diffusion/base_oxygen_diffusion.rst
@@ -0,0 +1,5 @@
+Base Model
+==========
+
+.. autoclass:: pybamm.oxygen_diffusion.BaseModel
+    :members:
diff --git a/docs/source/models/submodels/oxygen_diffusion/full_oxygen_diffusion.rst b/docs/source/models/submodels/oxygen_diffusion/full_oxygen_diffusion.rst
new file mode 100644
index 000000000..ad9e5e5a1
--- /dev/null
+++ b/docs/source/models/submodels/oxygen_diffusion/full_oxygen_diffusion.rst
@@ -0,0 +1,5 @@
+Full Model
+==========
+
+.. autoclass:: pybamm.oxygen_diffusion.Full
+    :members:
diff --git a/docs/source/models/submodels/oxygen_diffusion/index.rst b/docs/source/models/submodels/oxygen_diffusion/index.rst
new file mode 100644
index 000000000..8a2474751
--- /dev/null
+++ b/docs/source/models/submodels/oxygen_diffusion/index.rst
@@ -0,0 +1,9 @@
+Oxygen Diffusion
+================
+
+.. toctree::
+
+  base_oxygen_diffusion
+  no_oxygen
+  leading_oxygen_diffusion
+  full_oxygen_diffusion
diff --git a/docs/source/models/submodels/oxygen_diffusion/leading_oxygen_diffusion.rst b/docs/source/models/submodels/oxygen_diffusion/leading_oxygen_diffusion.rst
new file mode 100644
index 000000000..32b639786
--- /dev/null
+++ b/docs/source/models/submodels/oxygen_diffusion/leading_oxygen_diffusion.rst
@@ -0,0 +1,5 @@
+Leading Order Model
+===================
+
+.. autoclass:: pybamm.oxygen_diffusion.LeadingOrder
+    :members:
diff --git a/docs/source/models/submodels/oxygen_diffusion/no_oxygen.rst b/docs/source/models/submodels/oxygen_diffusion/no_oxygen.rst
new file mode 100644
index 000000000..ad938ec7c
--- /dev/null
+++ b/docs/source/models/submodels/oxygen_diffusion/no_oxygen.rst
@@ -0,0 +1,5 @@
+No Oxygen
+=========
+
+.. autoclass:: pybamm.oxygen_diffusion.NoOxygen
+    :members:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 7b167e7cb..a51de97bb 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -155,12 +155,13 @@ from .models.simple_ode_model import SimpleODEModel
 from .models.submodels.base_submodel import BaseSubModel
 
 from .models.submodels import (
+    convection,
+    current_collector,
     electrolyte,
     electrode,
-    particle,
-    current_collector,
-    convection,
     interface,
+    oxygen_diffusion,
+    particle,
     porosity,
     thermal,
 )
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index f196ddad5..67a3a1b4a 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -341,6 +341,8 @@ class Discretisation(object):
         :class:`pybamm.Matrix`
             The mass matrix
         """
+        pybamm.logger.debug("Create mass matrix")
+
         # Create list of mass matrices for each equation to be put into block
         # diagonal mass matrix for the model
         mass_list = []
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 873f99145..208122b89 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -348,15 +348,15 @@ class Multiplication(BinaryOperator):
 
         # anything multiplied by a scalar zero returns a scalar zero
         if is_scalar_zero(left):
-            if isinstance(right, pybamm.Array):
-                return pybamm.Array(np.zeros(right.shape))
-            else:
+            if right.shape_for_testing == ():
                 return pybamm.Scalar(0)
-        if is_scalar_zero(right):
-            if isinstance(left, pybamm.Array):
-                return pybamm.Array(np.zeros(left.shape))
             else:
+                return pybamm.Array(np.zeros(right.shape))
+        if is_scalar_zero(right):
+            if left.shape_for_testing == ():
                 return pybamm.Scalar(0)
+            else:
+                return pybamm.Array(np.zeros(left.shape))
 
         # if one of the children is a zero matrix, we have to be careful about shapes
         if is_matrix_zero(left) or is_matrix_zero(right):
@@ -477,11 +477,17 @@ class Division(BinaryOperator):
 
         # zero divided by anything returns zero
         if is_scalar_zero(left):
-            return pybamm.Scalar(0)
+            if right.shape_for_testing == ():
+                return pybamm.Scalar(0)
+            else:
+                return pybamm.Array(np.zeros(right.shape))
 
         # anything divided by zero returns inf
         if is_scalar_zero(right):
-            return pybamm.Scalar(np.inf)
+            if left.shape_for_testing == ():
+                return pybamm.Scalar(np.inf)
+            else:
+                return pybamm.Array(np.inf * np.ones(left.shape))
 
         # anything divided by one is itself
         if is_one(right):
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 9a4ce05c2..66497dfd6 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -530,6 +530,7 @@ class Symbol(anytree.NodeMixin):
             and not self.has_symbol_of_class(pybamm.Divergence)
             and not self.has_symbol_of_class(pybamm.IndefiniteIntegral)
             and not self.has_symbol_of_class(pybamm.Inner)
+            and not self.has_symbol_of_class(pybamm.Index)
         )
 
     def has_symbol_of_class(self, symbol_class):
@@ -585,8 +586,7 @@ class Symbol(anytree.NodeMixin):
         """
         Shape of an object for cases where it cannot be evaluated directly. If a symbol
         cannot be evaluated directly (e.g. it is a `Variable` or `Parameter`), it is
-        instead given an arbitrary domain-dependent shape from the dictionary
-        `pybamm.DOMAIN_SIZES_FOR_TESTING` (note that this only works for some domains)
+        instead given an arbitrary domain-dependent shape.
         """
         evaluated_self = self.evaluate_for_shape()
         if isinstance(evaluated_self, numbers.Number):
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index f41676113..86b9f3611 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -38,9 +38,8 @@ class BaseModel(object):
 
     """
 
-    def __init__(self):
-        # Default name
-        self.name = "Unnamed Model"
+    def __init__(self, name="Unnamed model"):
+        self.name = name
 
         # Initialise empty model
         self._rhs = {}
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 9437139f6..77c1f60a0 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -14,8 +14,8 @@ class BaseBatteryModel(pybamm.BaseModel):
     **Extends:** :class:`pybamm.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__()
+    def __init__(self, options=None, name="Unnamed battery model"):
+        super().__init__(name)
         self._extra_options = options
         self.set_standard_output_variables()
         self.submodels = OrderedDict()  # ordered dict not default in 3.5
@@ -290,27 +290,45 @@ class BaseBatteryModel(pybamm.BaseModel):
     def build_model(self):
 
         # Get the fundamental variables
-        for name, submodel in self.submodels.items():
-            pybamm.logger.debug("Getting fundamental variables for {}".format(name))
+        for submodel_name, submodel in self.submodels.items():
+            pybamm.logger.debug(
+                "Getting fundamental variables for {} ({})".format(
+                    submodel_name, self.name
+                )
+            )
             self.variables.update(submodel.get_fundamental_variables())
 
         # Get coupled variables
-        for name, submodel in self.submodels.items():
-            pybamm.logger.debug("Getting coupled variables for {}".format(name))
+        for submodel_name, submodel in self.submodels.items():
+            pybamm.logger.debug(
+                "Getting coupled variables for {} ({})".format(submodel_name, self.name)
+            )
             self.variables.update(submodel.get_coupled_variables(self.variables))
 
             # Set model equations
-        for name, submodel in self.submodels.items():
-            pybamm.logger.debug("Setting rhs for {}".format(name))
+        for submodel_name, submodel in self.submodels.items():
+            pybamm.logger.debug(
+                "Setting rhs for {} ({})".format(submodel_name, self.name)
+            )
             submodel.set_rhs(self.variables)
-            pybamm.logger.debug("Setting algebraic for {}".format(name))
+            pybamm.logger.debug(
+                "Setting algebraic for {} ({})".format(submodel_name, self.name)
+            )
             submodel.set_algebraic(self.variables)
-            pybamm.logger.debug("Setting boundary conditions for {}".format(name))
+            pybamm.logger.debug(
+                "Setting boundary conditions for {} ({})".format(
+                    submodel_name, self.name
+                )
+            )
             submodel.set_boundary_conditions(self.variables)
-            pybamm.logger.debug("Setting initial conditions for {}".format(name))
+            pybamm.logger.debug(
+                "Setting initial conditions for {} ({})".format(
+                    submodel_name, self.name
+                )
+            )
             submodel.set_initial_conditions(self.variables)
             submodel.set_events(self.variables)
-            pybamm.logger.debug("Updating {}".format(name))
+            pybamm.logger.debug("Updating {} ({})".format(submodel_name, self.name))
             self.update(submodel)
 
         pybamm.logger.debug("Setting voltage variables")
diff --git a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
index 8d9ded176..e6f77dfa7 100644
--- a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
+++ b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
@@ -16,8 +16,8 @@ class BaseModel(pybamm.BaseBatteryModel):
 
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
+    def __init__(self, options=None, name="Unnamed lead-acid model"):
+        super().__init__(options, name)
         self.param = pybamm.standard_parameters_lead_acid
 
     @property
@@ -81,3 +81,31 @@ class BaseModel(pybamm.BaseBatteryModel):
                 "Discharge capacity [A.h]": I * pybamm.t * time_scale / 3600,
             }
         )
+
+    def set_reactions(self):
+
+        # Should probably refactor as this is a bit clunky at the moment
+        # Maybe each reaction as a Reaction class so we can just list names of classes
+        param = self.param
+        icd = " interfacial current density"
+        self.reactions = {
+            "main": {
+                "Negative": {"s": param.s_n, "aj": "Negative electrode" + icd},
+                "Positive": {"s": param.s_p, "aj": "Positive electrode" + icd},
+            }
+        }
+        if "oxygen" in self.options["side reactions"]:
+            self.reactions["oxygen"] = {
+                "Negative": {
+                    "s": -(param.s_plus_Ox + param.t_plus),
+                    "s_ox": -param.s_ox_Ox,
+                    "aj": "Negative electrode oxygen" + icd,
+                },
+                "Positive": {
+                    "s": -(param.s_plus_Ox + param.t_plus),
+                    "s_ox": -param.s_ox_Ox,
+                    "aj": "Positive electrode oxygen" + icd,
+                },
+            }
+            self.reactions["main"]["Negative"]["s_ox"] = 0
+            self.reactions["main"]["Positive"]["s_ox"] = 0
diff --git a/pybamm/models/full_battery_models/lead_acid/composite.py b/pybamm/models/full_battery_models/lead_acid/composite.py
index e9791a402..e8ae80809 100644
--- a/pybamm/models/full_battery_models/lead_acid/composite.py
+++ b/pybamm/models/full_battery_models/lead_acid/composite.py
@@ -19,10 +19,10 @@ class Composite(BaseModel):
     **Extends:** :class:`pybamm.lead_acid.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
-        self.name = "Composite model"
+    def __init__(self, options=None, name="Composite model"):
+        super().__init__(options, name)
 
+        self.set_reactions()
         self.set_current_collector_submodel()
         self.set_interfacial_submodel()
         self.set_porosity_submodel()
@@ -80,7 +80,9 @@ class Composite(BaseModel):
 
         electrolyte = pybamm.electrolyte.stefan_maxwell
 
-        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(self.param)
+        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
+            self.param, self.reactions
+        )
 
         self.submodels[
             "electrolyte conductivity"
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index a3b6f05de..c488e6ca9 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -18,9 +18,8 @@ class LOQS(BaseModel):
     **Extends:** :class:`pybamm.lead_acid.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
-        self.name = "LOQS model"
+    def __init__(self, options=None, name="LOQS model"):
+        super().__init__(options, name)
         self.use_jacobian = False
 
         self.set_reactions()
@@ -32,26 +31,10 @@ class LOQS(BaseModel):
         self.set_electrolyte_submodel()
         self.set_positive_electrode_submodel()
         self.set_thermal_submodel()
+        self.set_side_reaction_submodels()
 
         self.build_model()
 
-    def set_reactions(self):
-
-        # Should probably refactor as this is a bit clunky at the moment
-        # Maybe each reaction as a Reaction class so we can just list names of classes
-        self.reactions = {
-            "main": {
-                "neg": {
-                    "s_plus": self.param.s_n,
-                    "j": "Average negative electrode interfacial current density",
-                },
-                "pos": {
-                    "s_plus": self.param.s_p,
-                    "j": "Average positive electrode interfacial current density",
-                },
-            }
-        }
-
     def set_current_collector_submodel(self):
 
         self.submodels["current collector"] = pybamm.current_collector.Uniform(
@@ -112,20 +95,43 @@ class LOQS(BaseModel):
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     domain.lower() + " electrolyte conductivity"
-                ] = surf_form.LeadingOrderDifferential(self.param, domain)
+                ] = surf_form.LeadingOrderDifferential(
+                    self.param, domain, self.reactions
+                )
         elif self.options["surface form"] == "algebraic":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     domain.lower() + " electrolyte conductivity"
-                ] = surf_form.LeadingOrderAlgebraic(self.param, domain)
-
-        else:
-            raise pybamm.OptionError("'capacitance' must be either 'True' or 'False'")
+                ] = surf_form.LeadingOrderAlgebraic(self.param, domain, self.reactions)
 
         self.submodels["electrolyte diffusion"] = electrolyte.diffusion.LeadingOrder(
             self.param, self.reactions
         )
 
+    def set_side_reaction_submodels(self):
+        if "oxygen" in self.options["side reactions"]:
+            self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.LeadingOrder(
+                self.param, self.reactions
+            )
+            self.submodels[
+                "positive oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.ForwardTafel(self.param, "Positive")
+            self.submodels[
+                "negative oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.LeadingOrderDiffusionLimited(
+                self.param, "Negative"
+            )
+        else:
+            self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.NoOxygen(
+                self.param
+            )
+            self.submodels[
+                "positive oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.NoReaction(self.param, "Positive")
+            self.submodels[
+                "negative oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.NoReaction(self.param, "Negative")
+
     @property
     def default_spatial_methods(self):
         # ODEs only in the macroscale, so use base spatial method
diff --git a/pybamm/models/full_battery_models/lead_acid/newman_tiedemann.py b/pybamm/models/full_battery_models/lead_acid/newman_tiedemann.py
index 36f942cbd..2a90f3e4b 100644
--- a/pybamm/models/full_battery_models/lead_acid/newman_tiedemann.py
+++ b/pybamm/models/full_battery_models/lead_acid/newman_tiedemann.py
@@ -18,10 +18,10 @@ class NewmanTiedemann(BaseModel):
     **Extends:** :class:`pybamm.lead_acid.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
-        self.name = "Newman-Tiedemann model"
+    def __init__(self, options=None, name="Newman-Tiedemann model"):
+        super().__init__(options, name)
 
+        self.set_reactions()
         self.set_current_collector_submodel()
         self.set_interfacial_submodel()
         self.set_porosity_submodel()
@@ -30,6 +30,7 @@ class NewmanTiedemann(BaseModel):
         self.set_negative_electrode_submodel()
         self.set_positive_electrode_submodel()
         self.set_thermal_submodel()
+        self.set_side_reaction_submodels()
 
         self.build_model()
 
@@ -57,14 +58,14 @@ class NewmanTiedemann(BaseModel):
 
     def set_negative_electrode_submodel(self):
         if self.options["surface form"] is False:
-            submodel = pybamm.electrode.ohm.Full(self.param, "Negative")
+            submodel = pybamm.electrode.ohm.Full(self.param, "Negative", self.reactions)
         else:
             submodel = pybamm.electrode.ohm.SurfaceForm(self.param, "Negative")
         self.submodels["negative electrode"] = submodel
 
     def set_positive_electrode_submodel(self):
         if self.options["surface form"] is False:
-            submodel = pybamm.electrode.ohm.Full(self.param, "Positive")
+            submodel = pybamm.electrode.ohm.Full(self.param, "Positive", self.reactions)
         else:
             submodel = pybamm.electrode.ohm.SurfaceForm(self.param, "Positive")
         self.submodels["positive electrode"] = submodel
@@ -74,22 +75,48 @@ class NewmanTiedemann(BaseModel):
         electrolyte = pybamm.electrolyte.stefan_maxwell
         surf_form = electrolyte.conductivity.surface_potential_form
 
-        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(self.param)
+        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
+            self.param, self.reactions
+        )
 
         if self.options["surface form"] is False:
             self.submodels["electrolyte conductivity"] = electrolyte.conductivity.Full(
-                self.param
+                self.param, self.reactions
             )
         elif self.options["surface form"] == "differential":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     domain.lower() + " electrolyte conductivity"
-                ] = surf_form.FullDifferential(self.param, domain)
+                ] = surf_form.FullDifferential(self.param, domain, self.reactions)
         elif self.options["surface form"] == "algebraic":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     domain.lower() + " electrolyte conductivity"
-                ] = surf_form.FullAlgebraic(self.param, domain)
+                ] = surf_form.FullAlgebraic(self.param, domain, self.reactions)
+
+    def set_side_reaction_submodels(self):
+        if "oxygen" in self.options["side reactions"]:
+            self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.Full(
+                self.param, self.reactions
+            )
+            self.submodels[
+                "positive oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.ForwardTafel(self.param, "Positive")
+            self.submodels[
+                "negative oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.FullDiffusionLimited(
+                self.param, "Negative"
+            )
+        else:
+            self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.NoOxygen(
+                self.param
+            )
+            self.submodels[
+                "positive oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.NoReaction(self.param, "Positive")
+            self.submodels[
+                "negative oxygen interface"
+            ] = pybamm.interface.lead_acid_oxygen.NoReaction(self.param, "Negative")
 
     @property
     def default_solver(self):
diff --git a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index 16e7a27d2..dd4565da6 100644
--- a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -13,8 +13,8 @@ class BaseModel(pybamm.BaseBatteryModel):
 
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
+    def __init__(self, options=None, name="Unnamed lithium-ion model"):
+        super().__init__(options, name)
         self.param = pybamm.standard_parameters_lithium_ion
 
     def set_standard_output_variables(self):
@@ -67,3 +67,21 @@ class BaseModel(pybamm.BaseBatteryModel):
                 "r_p [m]": var.r_p * param.R_p,
             }
         )
+
+    def set_reactions(self):
+
+        # Should probably refactor as this is a bit clunky at the moment
+        # Maybe each reaction as a Reaction class so we can just list names of classes
+        icd = " interfacial current density"
+        self.reactions = {
+            "main": {
+                "Negative": {
+                    "s": 1 - self.param.t_plus,
+                    "aj": "Negative electrode" + icd,
+                },
+                "Positive": {
+                    "s": 1 - self.param.t_plus,
+                    "aj": "Positive electrode" + icd,
+                },
+            }
+        }
diff --git a/pybamm/models/full_battery_models/lithium_ion/dfn.py b/pybamm/models/full_battery_models/lithium_ion/dfn.py
index d00be5bb4..4a2c801fb 100644
--- a/pybamm/models/full_battery_models/lithium_ion/dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/dfn.py
@@ -12,10 +12,10 @@ class DFN(BaseModel):
     **Extends:** :class:`pybamm.lithium_ion.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
-        self.name = "Doyle-Fuller-Newman model"
+    def __init__(self, options=None, name="Doyle-Fuller-Newman model"):
+        super().__init__(options, name)
 
+        self.set_reactions()
         self.set_current_collector_submodel()
         self.set_porosity_submodel()
         self.set_convection_submodel()
@@ -62,10 +62,10 @@ class DFN(BaseModel):
     def set_solid_submodel(self):
 
         self.submodels["negative electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Negative"
+            self.param, "Negative", self.reactions
         )
         self.submodels["positive electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Positive"
+            self.param, "Positive", self.reactions
         )
 
     def set_electrolyte_submodel(self):
@@ -73,9 +73,11 @@ class DFN(BaseModel):
         electrolyte = pybamm.electrolyte.stefan_maxwell
 
         self.submodels["electrolyte conductivity"] = electrolyte.conductivity.Full(
-            self.param
+            self.param, self.reactions
+        )
+        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
+            self.param, self.reactions
         )
-        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(self.param)
 
     @property
     def default_geometry(self):
diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py
index 0c499be8c..81300045d 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spm.py
@@ -11,9 +11,8 @@ class SPM(BaseModel):
     **Extends:** :class:`pybamm.lithium_ion.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
-        self.name = "Single particle model"
+    def __init__(self, options=None, name="Single Particle Model"):
+        super().__init__(options, name)
 
         self.set_current_collector_submodel()
         self.set_porosity_submodel()
diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py
index ea8738f02..fbf6efb0f 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spme.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spme.py
@@ -11,10 +11,10 @@ class SPMe(BaseModel):
     **Extends:** :class:`pybamm.lithium_ion.BaseModel`
     """
 
-    def __init__(self, options=None):
-        super().__init__(options)
-        self.name = "Single Particle Model with electrolyte"
+    def __init__(self, options=None, name="Single Particle Model with electrolyte"):
+        super().__init__(options, name)
 
+        self.set_reactions()
         self.set_current_collector_submodel()
         self.set_porosity_submodel()
         self.set_convection_submodel()
@@ -103,7 +103,9 @@ class SPMe(BaseModel):
         self.submodels[
             "electrolyte conductivity"
         ] = electrolyte.conductivity.CombinedOrder(self.param)
-        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(self.param)
+        self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
+            self.param, self.reactions
+        )
 
     @property
     def default_geometry(self):
diff --git a/pybamm/models/reaction_diffusion.py b/pybamm/models/reaction_diffusion.py
index fbd2ac520..de7204c12 100644
--- a/pybamm/models/reaction_diffusion.py
+++ b/pybamm/models/reaction_diffusion.py
@@ -12,15 +12,7 @@ class ReactionDiffusionModel(pybamm.BaseBatteryModel):
     """
 
     def __init__(self, options=None):
-
-        if not options:
-            options = {}
-
-        options.update(
-            {"Voltage": "Off"}
-        )  # annoying option only for reaction diffusion
-        super().__init__(options)
-        self.name = "Reaction diffusion model"
+        super().__init__(options, "Reaction diffusion model")
         # NOTE: set_standard_output_variables sets the lead acid timescales,
         # so if paramaters are changed here the timescale in the method
         # set_standard_output_variables may need to be altered
@@ -31,6 +23,7 @@ class ReactionDiffusionModel(pybamm.BaseBatteryModel):
         self.param.epsilon_s = pybamm.Scalar(1)
         self.param.epsilon_p = pybamm.Scalar(1)
 
+        self.set_reactions()
         self.set_current_collector_submodel()
         self.set_convection_submodel()
         self.set_porosity_submodel()
@@ -54,7 +47,7 @@ class ReactionDiffusionModel(pybamm.BaseBatteryModel):
     def set_electrolyte_submodel(self):
         electrolyte = pybamm.electrolyte.stefan_maxwell
         self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
-            self.param
+            self.param, self.reactions
         )
 
     def set_interfacial_submodel(self):
@@ -99,3 +92,16 @@ class ReactionDiffusionModel(pybamm.BaseBatteryModel):
                 "Discharge capacity [A.h]": I * pybamm.t * time_scale / 3600,
             }
         )
+
+    def set_reactions(self):
+
+        # Should probably refactor as this is a bit clunky at the moment
+        # Maybe each reaction as a Reaction class so we can just list names of classes
+        param = self.param
+        icd = " interfacial current density"
+        self.reactions = {
+            "main": {
+                "Negative": {"s": param.s_n, "aj": "Negative electrode" + icd},
+                "Positive": {"s": param.s_p, "aj": "Positive electrode" + icd},
+            }
+        }
diff --git a/pybamm/models/simple_ode_model.py b/pybamm/models/simple_ode_model.py
index 9084133ad..705ecbda4 100644
--- a/pybamm/models/simple_ode_model.py
+++ b/pybamm/models/simple_ode_model.py
@@ -14,8 +14,7 @@ class SimpleODEModel(pybamm.BaseBatteryModel):
     """
 
     def __init__(self):
-        super().__init__()
-        self.name = "Simple ODE Model"
+        super().__init__(None, "Simple ODE Model")
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         # Create variables: domain is explicitly empty since these variables are only
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index cd3b1ccdf..26f41a2af 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -46,7 +46,7 @@ class BaseSubModel:
         symbols.
     """
 
-    def __init__(self, param, domain=None):
+    def __init__(self, param, domain=None, reactions=None):
         super().__init__()
         self.param = param
         # Initialise empty variables (to avoid overwriting with 'None')
@@ -60,6 +60,7 @@ class BaseSubModel:
 
         self.domain = domain
         self.set_domain_for_broadcast()
+        self.reactions = reactions
 
     @property
     def domain(self):
diff --git a/pybamm/models/submodels/electrode/base_electrode.py b/pybamm/models/submodels/electrode/base_electrode.py
index 3afccfe88..1fb3c1b67 100644
--- a/pybamm/models/submodels/electrode/base_electrode.py
+++ b/pybamm/models/submodels/electrode/base_electrode.py
@@ -17,8 +17,8 @@ class BaseElectrode(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions=None):
+        super().__init__(param, domain, reactions)
 
     def _get_standard_potential_variables(self, phi_s):
         """
diff --git a/pybamm/models/submodels/electrode/ohm/base_ohm.py b/pybamm/models/submodels/electrode/ohm/base_ohm.py
index 6b4fdef9d..f7755d550 100644
--- a/pybamm/models/submodels/electrode/ohm/base_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/base_ohm.py
@@ -20,8 +20,8 @@ class BaseModel(BaseElectrode):
     **Extends:** :class:`pybamm.electrode.BaseElectrode`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions=None):
+        super().__init__(param, domain, reactions)
 
     @property
     def default_solver(self):
diff --git a/pybamm/models/submodels/electrode/ohm/full_ohm.py b/pybamm/models/submodels/electrode/ohm/full_ohm.py
index f7625bf1e..93d34c0e5 100644
--- a/pybamm/models/submodels/electrode/ohm/full_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/full_ohm.py
@@ -19,8 +19,8 @@ class Full(BaseModel):
     **Extends:** :class:`pybamm.electrode.ohm.BaseModel`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def get_fundamental_variables(self):
 
@@ -59,9 +59,12 @@ class Full(BaseModel):
 
         phi_s = variables[self.domain + " electrode potential"]
         i_s = variables[self.domain + " electrode current density"]
-        j = variables[self.domain + " electrode interfacial current density"]
+        sum_j = sum(
+            variables[reaction[self.domain]["aj"]]
+            for reaction in self.reactions.values()
+        )
 
-        self.algebraic[phi_s] = pybamm.div(i_s) + j
+        self.algebraic[phi_s] = pybamm.div(i_s) + sum_j
 
     def set_boundary_conditions(self, variables):
 
diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py
index c34513926..d8ce9927f 100644
--- a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py
@@ -17,8 +17,9 @@ class BaseElectrolyteConductivity(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, domain=None):
+    def __init__(self, param, domain=None, reactions=None):
         super().__init__(param, domain)
+        self.reactions = reactions
 
     def _get_standard_potential_variables(self, phi_e, phi_e_av):
         """
diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
index 365bed078..6c2b57696 100644
--- a/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
@@ -16,8 +16,8 @@ class BaseElectrolyteDiffusion(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, reactions=None):
+        super().__init__(param, reactions=reactions)
 
     def _get_standard_concentration_variables(self, c_e):
         """
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py
index 7f6ff1975..198d6bf1c 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py
@@ -17,5 +17,5 @@ class BaseModel(BaseElectrolyteConductivity):
     **Extends:** :class:`pybamm.electrolyte.BaseElectrolyteConductivity`
     """
 
-    def __init__(self, param, domain=None):
-        super().__init__(param, domain)
+    def __init__(self, param, domain=None, reactions=None):
+        super().__init__(param, domain, reactions)
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
index a5b9de666..d4d8c34eb 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
@@ -20,8 +20,8 @@ class Full(BaseModel):
     **Extends:** :class:`pybamm.BaseStefanMaxwellConductivity`
     """
 
-    def __init__(self, param, domain=None):
-        super().__init__(param, domain)
+    def __init__(self, param, reactions):
+        super().__init__(param, reactions=reactions)
 
     def get_fundamental_variables(self):
         phi_e = pybamm.standard_variables.phi_e
@@ -47,9 +47,16 @@ class Full(BaseModel):
     def set_algebraic(self, variables):
         phi_e = variables["Electrolyte potential"]
         i_e = variables["Electrolyte current density"]
-        j = variables["Interfacial current density"]
+        sum_j = sum(
+            pybamm.Concatenation(
+                variables[reaction["Negative"]["aj"]],
+                pybamm.Broadcast(0, "separator"),
+                variables[reaction["Positive"]["aj"]],
+            )
+            for reaction in self.reactions.values()
+        )
 
-        self.algebraic = {phi_e: pybamm.div(i_e) - j}
+        self.algebraic = {phi_e: pybamm.div(i_e) - sum_j}
 
     def set_boundary_conditions(self, variables):
         phi_e = variables["Electrolyte potential"]
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/leading_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/leading_stefan_maxwell_conductivity.py
index 45fd8b089..3ec8954f2 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/leading_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/leading_stefan_maxwell_conductivity.py
@@ -19,8 +19,8 @@ class LeadingOrder(BaseModel):
     **Extends:** :class:`pybamm.BaseStefanMaxwellConductivity`
     """
 
-    def __init__(self, param, domain=None):
-        super().__init__(param, domain)
+    def __init__(self, param, domain=None, reactions=None):
+        super().__init__(param, domain, reactions)
 
     def get_coupled_variables(self, variables):
         ocp_n_av = variables["Average negative electrode open circuit potential"]
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
index 39ac1661f..a19613c07 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
@@ -21,8 +21,8 @@ class BaseModel(BaseStefanMaxwellConductivity):
     **Extends:** :class:`pybamm.electrolyte.stefan_maxwell.conductivity.BaseModel`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def get_fundamental_variables(self):
         if self.domain == "Negative":
@@ -31,8 +31,6 @@ class BaseModel(BaseStefanMaxwellConductivity):
             return {}
         elif self.domain == "Positive":
             delta_phi = pybamm.standard_variables.delta_phi_p
-        else:
-            raise pybamm.DomainError
 
         variables = self._get_standard_surface_potential_difference_variables(delta_phi)
 
@@ -47,8 +45,6 @@ class BaseModel(BaseStefanMaxwellConductivity):
             delta_phi_e_init = self.param.U_n(self.param.c_n_init)
         elif self.domain == "Positive":
             delta_phi_e_init = self.param.U_p(self.param.c_p_init)
-        else:
-            raise pybamm.DomainError
 
         self.initial_conditions = {delta_phi_e: delta_phi_e_init}
 
@@ -103,9 +99,6 @@ class BaseModel(BaseStefanMaxwellConductivity):
             lbc_c_e = (c_e_flux, "Neumann")
             rbc_c_e = (pybamm.Scalar(0), "Neumann")
 
-        else:
-            raise pybamm.DomainError
-
         # TODO: check if we still need the boundary conditions for c_e, once we have
         # internal boundary conditions
         self.boundary_conditions = {
@@ -241,8 +234,8 @@ class FullAlgebraic(BaseModel):
      **Extends:** :class:`pybamm.electrolyte.stefan_maxwell.conductivity.surface_potential_form.BaseFull`
     """  # noqa: E501
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def set_algebraic(self, variables):
         if self.domain == "Separator":
@@ -250,9 +243,11 @@ class FullAlgebraic(BaseModel):
 
         delta_phi = variables[self.domain + " electrode surface potential difference"]
         i_e = variables[self.domain + " electrolyte current density"]
-        j = variables[self.domain + " electrode interfacial current density"]
-
-        self.algebraic[delta_phi] = pybamm.div(i_e) - j
+        sum_j = sum(
+            variables[reaction[self.domain]["aj"]]
+            for reaction in self.reactions.values()
+        )
+        self.algebraic[delta_phi] = pybamm.div(i_e) - sum_j
 
 
 class FullDifferential(BaseModel):
@@ -270,8 +265,8 @@ class FullDifferential(BaseModel):
 
     """  # noqa: E501
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def set_rhs(self, variables):
         if self.domain == "Separator":
@@ -284,6 +279,9 @@ class FullDifferential(BaseModel):
 
         delta_phi = variables[self.domain + " electrode surface potential difference"]
         i_e = variables[self.domain + " electrolyte current density"]
-        j = variables[self.domain + " electrode interfacial current density"]
+        sum_j = sum(
+            variables[reaction[self.domain]["aj"]]
+            for reaction in self.reactions.values()
+        )
 
-        self.rhs[delta_phi] = 1 / C_dl * (pybamm.div(i_e) - j)
+        self.rhs[delta_phi] = 1 / C_dl * (pybamm.div(i_e) - sum_j)
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py
index a07163fdc..c1260cddd 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py
@@ -19,8 +19,8 @@ class BaseLeadingOrderSurfaceForm(LeadingOrder):
     **Extends:** :class:`pybamm.electrolyte.stefan_maxwell.conductivity.surface_potential_form.BaseModel`
     """  # noqa: E501
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def get_fundamental_variables(self):
 
@@ -30,8 +30,6 @@ class BaseLeadingOrderSurfaceForm(LeadingOrder):
             return {}
         elif self.domain == "Positive":
             delta_phi = pybamm.standard_variables.delta_phi_p_av
-        else:
-            raise pybamm.DomainError
 
         variables = self._get_standard_surface_potential_difference_variables(delta_phi)
         return variables
@@ -49,9 +47,6 @@ class BaseLeadingOrderSurfaceForm(LeadingOrder):
         elif self.domain == "Positive":
             delta_phi_init = self.param.U_p(self.param.c_p_init)
 
-        else:
-            raise pybamm.DomainError
-
         self.initial_conditions = {delta_phi: delta_phi_init}
 
     def get_coupled_variables(self, variables):
@@ -61,10 +56,7 @@ class BaseLeadingOrderSurfaceForm(LeadingOrder):
             delta_phi_n_av = variables[
                 "Average negative electrode surface potential difference"
             ]
-
-            # Potential
             phi_e_av = -delta_phi_n_av
-
             return self._get_coupled_variables_from_potential(variables, phi_e_av)
 
         else:
@@ -86,8 +78,8 @@ class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
 
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def set_rhs(self, variables):
         if self.domain == "Separator":
@@ -95,8 +87,12 @@ class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
 
         param = self.param
 
-        j = variables[self.domain + " electrode interfacial current density"].orphans[0]
-        j_tot_av = variables[
+        sum_j = sum(
+            variables[reaction[self.domain]["aj"]].orphans[0]
+            for reaction in self.reactions.values()
+        )
+
+        sum_j_av = variables[
             "Average "
             + self.domain.lower()
             + " electrode total interfacial current density"
@@ -110,7 +106,7 @@ class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
         elif self.domain == "Positive":
             C_dl = param.C_dl_p
 
-        self.rhs[delta_phi] = 1 / C_dl * (j_tot_av - j)
+        self.rhs[delta_phi] = 1 / C_dl * (sum_j_av - sum_j)
 
 
 class LeadingOrderAlgebraic(BaseLeadingOrderSurfaceForm):
@@ -127,15 +123,19 @@ class LeadingOrderAlgebraic(BaseLeadingOrderSurfaceForm):
     **Extends:** :class:`BaseLeadingOrderSurfaceForm`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, reactions):
+        super().__init__(param, domain, reactions)
 
     def set_algebraic(self, variables):
         if self.domain == "Separator":
             return
 
-        j = variables[self.domain + " electrode interfacial current density"].orphans[0]
-        j_tot_av = variables[
+        sum_j = sum(
+            variables[reaction[self.domain]["aj"]].orphans[0]
+            for reaction in self.reactions.values()
+        )
+
+        sum_j_av = variables[
             "Average "
             + self.domain.lower()
             + " electrode total interfacial current density"
@@ -144,4 +144,4 @@ class LeadingOrderAlgebraic(BaseLeadingOrderSurfaceForm):
             "Average " + self.domain.lower() + " electrode surface potential difference"
         ]
 
-        self.algebraic[delta_phi] = j_tot_av - j
+        self.algebraic[delta_phi] = sum_j_av - sum_j
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py
index fc7354131..0bb8ec46e 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py
@@ -17,5 +17,5 @@ class BaseModel(BaseElectrolyteDiffusion):
     **Extends:** :class:`pybamm.electrolyte.BaseElectrolyteDiffusion`
     """
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, reactions=None):
+        super().__init__(param, reactions)
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py
index f6fdb14af..fd13944bf 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py
@@ -7,9 +7,7 @@ from .base_stefan_maxwell_diffusion import BaseModel
 
 
 class ConstantConcentration(BaseModel):
-    """Class for conservation of mass in the electrolyte employing the
-    Stefan-Maxwell constitutive equations. (Leading refers to leading order
-    of asymptotic reduction)
+    """Class for constant concentration of electrolyte
 
     Parameters
     ----------
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
index 61b93fd16..31eddfef8 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
@@ -20,8 +20,8 @@ class Full(BaseModel):
     **Extends:** :class:`pybamm.electrolyte.stefan_maxwell.diffusion.BaseModel`
     """
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, reactions):
+        super().__init__(param, reactions)
 
     def get_fundamental_variables(self):
         c_e = pybamm.standard_variables.c_e
@@ -58,16 +58,23 @@ class Full(BaseModel):
         c_e = variables["Electrolyte concentration"]
         N_e = variables["Electrolyte flux"]
         # i_e = variables["Electrolyte current density"]
-        j = variables["Interfacial current density"]
 
         # TODO: check lead acid version in new form
-        source_term = param.s / param.gamma_e * j
         # source_term = ((param.s - param.t_plus) / param.gamma_e) * pybamm.div(i_e)
         # source_term = pybamm.div(i_e) / param.gamma_e  # lithium-ion
+        source_terms = sum(
+            pybamm.Concatenation(
+                reaction["Negative"]["s"] * variables[reaction["Negative"]["aj"]],
+                pybamm.Broadcast(0, "separator"),
+                reaction["Positive"]["s"] * variables[reaction["Positive"]["aj"]],
+            )
+            / param.gamma_e
+            for reaction in self.reactions.values()
+        )
 
         self.rhs = {
             c_e: (1 / eps)
-            * (-pybamm.div(N_e) / param.C_e + source_term - c_e * deps_dt)
+            * (-pybamm.div(N_e) / param.C_e + source_terms - c_e * deps_dt)
         }
 
     def set_boundary_conditions(self, variables):
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py
index 00b5edd11..35d69a4b2 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py
@@ -15,14 +15,14 @@ class LeadingOrder(BaseModel):
     ----------
     param : parameter class
         The parameters to use for this submodel
-
+    reactions : dict
+        Dictionary of reaction terms
 
     **Extends:** :class:`pybamm.electrolyte.stefan_maxwell.diffusion.BaseModel`
     """
 
     def __init__(self, param, reactions):
-        super().__init__(param)
-        self.reactions = reactions
+        super().__init__(param, reactions)
 
     def get_fundamental_variables(self):
         c_e_av = pybamm.standard_variables.c_e_av
@@ -56,21 +56,16 @@ class LeadingOrder(BaseModel):
         deps_n_dt_av = variables["Average negative electrode porosity change"]
         deps_p_dt_av = variables["Average positive electrode porosity change"]
 
-        # TODO: ask tino about this bit
-        # the "j" component of reactions is now just a string referring to the variable
-        # so that the expression doesn't need to be set twice
         source_terms = sum(
-            param.l_n * rxn["neg"]["s_plus"] * variables[rxn["neg"]["j"]]
-            + param.l_p * rxn["pos"]["s_plus"] * variables[rxn["pos"]["j"]]
+            param.l_n
+            * rxn["Negative"]["s"]
+            * variables[rxn["Negative"]["aj"]].orphans[0]
+            + param.l_p
+            * rxn["Positive"]["s"]
+            * variables[rxn["Positive"]["aj"]].orphans[0]
             for rxn in self.reactions.values()
         )
 
-        # source_terms = sum(
-        #     param.l_n * rxn["neg"]["s_plus"] * rxn["neg"]["aj"]
-        #     + param.l_p * rxn["pos"]["s_plus"] * rxn["pos"]["aj"]
-        #     for rxn in self.reactions.values()
-        # )
-
         self.rhs = {
             c_e_av: 1
             / (param.l_n * eps_n_av + param.l_s * eps_s_av + param.l_p * eps_p_av)
diff --git a/pybamm/models/submodels/interface/__init__.py b/pybamm/models/submodels/interface/__init__.py
index 32ea1b15b..179aa1eb7 100644
--- a/pybamm/models/submodels/interface/__init__.py
+++ b/pybamm/models/submodels/interface/__init__.py
@@ -1,4 +1,5 @@
 from .base_interface import BaseInterface
 from .kinetics import BaseModel
 from . import lead_acid
+from . import lead_acid_oxygen
 from . import lithium_ion
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index 7b6ad6bad..5f4cd358c 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -21,6 +21,16 @@ class BaseInterface(pybamm.BaseSubModel):
     def __init__(self, param, domain):
         super().__init__(param, domain)
 
+    def _get_delta_phi_s(self, variables):
+        "Calculate delta_phi_s, and derived variables, using phi_s and phi_e"
+        phi_s = variables[self.domain + " electrode potential"]
+        phi_e = variables[self.domain + " electrolyte potential"]
+        delta_phi_s = phi_s - phi_e
+        variables.update(
+            self._get_standard_surface_potential_difference_variables(delta_phi_s)
+        )
+        return variables
+
     def _get_average_total_interfacial_current_density(self, variables):
         """
         Method to obtain the average total interfacial current density.
@@ -48,14 +58,23 @@ class BaseInterface(pybamm.BaseSubModel):
             j = pybamm.Broadcast(j, self.domain_for_broadcast, broadcast_type="primary")
 
         variables = {
-            self.domain + " electrode interfacial current density": j,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " interfacial current density": j,
             "Average "
             + self.domain.lower()
-            + " electrode interfacial current density": j_av,
-            self.domain + " interfacial current density [A.m-2]": i_typ * j,
+            + " electrode"
+            + self.reaction_name
+            + " interfacial current density": j_av,
+            self.domain
+            + self.reaction_name
+            + " interfacial current density [A.m-2]": i_typ * j,
             "Average "
             + self.domain.lower()
-            + " electrode interfacial current density [A.m-2]": i_typ * j_av,
+            + " electrode"
+            + self.reaction_name
+            + " interfacial current density [A.m-2]": i_typ * j_av,
         }
 
         return variables
@@ -84,12 +103,17 @@ class BaseInterface(pybamm.BaseSubModel):
         j_p = variables["Positive electrode interfacial current density"]
         j = pybamm.Concatenation(j_n, j_s, j_p)
 
-        variables.update(
-            {
+        if self.reaction_name == "":
+            variables = {
                 "Interfacial current density": j,
                 "Interfacial current density [A.m-2]": i_typ * j,
             }
-        )
+        else:
+            reaction_name = self.reaction_name[1:].capitalize()
+            variables = {
+                reaction_name + " interfacial current density": j,
+                reaction_name + " interfacial current density [A.m-2]": i_typ * j,
+            }
 
         return variables
 
@@ -106,14 +130,24 @@ class BaseInterface(pybamm.BaseSubModel):
             )
 
         variables = {
-            self.domain + " electrode exchange current density": j0,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " exchange current density": j0,
             "Average "
             + self.domain.lower()
-            + " electrode exchange current density": j0_av,
-            self.domain + " electrode exchange current density [A.m-2]": i_typ * j0,
+            + " electrode"
+            + self.reaction_name
+            + " exchange current density": j0_av,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " exchange current density [A.m-2]": i_typ * j0,
             "Average "
             + self.domain.lower()
-            + " electrode exchange current density [A.m-2]": i_typ * j0_av,
+            + " electrode"
+            + self.reaction_name
+            + " exchange current density [A.m-2]": i_typ * j0_av,
         }
 
         return variables
@@ -127,12 +161,17 @@ class BaseInterface(pybamm.BaseSubModel):
         j0_p = variables["Positive electrode exchange current density"]
         j0 = pybamm.Concatenation(j0_n, j0_s, j0_p)
 
-        variables.update(
-            {
+        if self.reaction_name == "":
+            variables = {
                 "Exchange current density": j0,
                 "Exchange current density [A.m-2]": i_typ * j0,
             }
-        )
+        else:
+            reaction_name = self.reaction_name[1:].capitalize()
+            variables = {
+                reaction_name + " exchange current density": j0,
+                reaction_name + " exchange current density [A.m-2]": i_typ * j0,
+            }
 
         return variables
 
@@ -149,14 +188,24 @@ class BaseInterface(pybamm.BaseSubModel):
             )
 
         variables = {
-            self.domain + " electrode reaction overpotential": eta_r,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " reaction overpotential": eta_r,
             "Average "
             + self.domain.lower()
-            + " electrode reaction overpotential": eta_r_av,
-            self.domain + " electrode reaction overpotential [V]": eta_r * pot_scale,
+            + " electrode"
+            + self.reaction_name
+            + " reaction overpotential": eta_r_av,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " reaction overpotential [V]": eta_r * pot_scale,
             "Average "
             + self.domain.lower()
-            + " electrode reaction overpotential [V]": eta_r_av * pot_scale,
+            + " electrode"
+            + self.reaction_name
+            + " reaction overpotential [V]": eta_r_av * pot_scale,
         }
 
         return variables
@@ -231,14 +280,24 @@ class BaseInterface(pybamm.BaseSubModel):
             ocp_av_dim = self.param.U_p_ref + self.param.potential_scale * ocp_av
 
         variables = {
-            self.domain + " electrode open circuit potential": ocp,
-            self.domain + " electrode open circuit potential [V]": ocp_dim,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " open circuit potential": ocp,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " open circuit potential [V]": ocp_dim,
             "Average "
             + self.domain.lower()
-            + " electrode open circuit potential": ocp_av,
+            + " electrode"
+            + self.reaction_name
+            + " open circuit potential": ocp_av,
             "Average "
             + self.domain.lower()
-            + " electrode open circuit potential [V]": ocp_av_dim,
+            + " electrode"
+            + self.reaction_name
+            + " open circuit potential [V]": ocp_av_dim,
             self.domain + " electrode entropic change": dUdT,
             "Average " + self.domain.lower() + " electrode entropic change": dUdT_av,
         }
diff --git a/pybamm/models/submodels/interface/diffusion_limited/__init__.py b/pybamm/models/submodels/interface/diffusion_limited/__init__.py
new file mode 100644
index 000000000..d6fc1404b
--- /dev/null
+++ b/pybamm/models/submodels/interface/diffusion_limited/__init__.py
@@ -0,0 +1,3 @@
+from .base_diffusion_limited import BaseModel
+from .leading_diffusion_limited import BaseLeadingOrderDiffusionLimited
+from .full_diffusion_limited import BaseFullDiffusionLimited
diff --git a/pybamm/models/submodels/interface/diffusion_limited/base_diffusion_limited.py b/pybamm/models/submodels/interface/diffusion_limited/base_diffusion_limited.py
new file mode 100644
index 000000000..e2fb02ab7
--- /dev/null
+++ b/pybamm/models/submodels/interface/diffusion_limited/base_diffusion_limited.py
@@ -0,0 +1,78 @@
+#
+# Leading-order diffusion limited kinetics
+#
+
+import pybamm
+from ..base_interface import BaseInterface
+
+
+class BaseModel(BaseInterface):
+    """
+    Leading-order submodel for diffusion-limited kinetics
+
+    Parameters
+    ----------
+    param :
+        model parameters
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
+
+
+    **Extends:** :class:`pybamm.interface.BaseInterface`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def get_coupled_variables(self, variables):
+        # Calculate delta_phi_s from phi_s and phi_e if it isn't already known
+        if self.domain + " electrode surface potential difference" not in variables:
+            variables = self._get_delta_phi_s(variables)
+        delta_phi_s = variables[self.domain + " electrode surface potential difference"]
+        # If delta_phi_s was broadcast, take only the orphan
+        if isinstance(delta_phi_s, pybamm.Broadcast):
+            delta_phi_s = delta_phi_s.orphans[0]
+
+        # Get exchange-current density
+        j0 = self._get_exchange_current_density(variables)
+        # Get open-circuit potential variables and reaction overpotential
+        ocp, dUdT = self._get_open_circuit_potential(variables)
+        eta_r = delta_phi_s - ocp
+
+        # Get interfacial current densities
+        j = self._get_diffusion_limited_current_density(variables)
+        j_tot_av = self._get_average_total_interfacial_current_density(variables)
+
+        variables.update(self._get_standard_interfacial_current_variables(j))
+        variables.update(
+            self._get_standard_total_interfacial_current_variables(j_tot_av)
+        )
+        variables.update(self._get_standard_exchange_current_variables(j0))
+        variables.update(self._get_standard_overpotential_variables(eta_r))
+        variables.update(self._get_standard_ocp_variables(ocp, dUdT))
+
+        if (
+            "Negative electrode" + self.reaction_name + " interfacial current density"
+            in variables
+            and "Positive electrode"
+            + self.reaction_name
+            + " interfacial current density"
+            in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+            variables.update(
+                self._get_standard_whole_cell_exchange_current_variables(variables)
+            )
+
+        return variables
+
+    def _get_exchange_current_density(self, variables):
+        raise NotImplementedError
+
+    def _get_open_circuit_potential(self, variables):
+        raise NotImplementedError
+
+    def _get_diffusion_limited_current_density(self, variables):
+        raise NotImplementedError
diff --git a/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py b/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py
new file mode 100644
index 000000000..88489e558
--- /dev/null
+++ b/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py
@@ -0,0 +1,33 @@
+#
+# Full diffusion limited kinetics
+#
+from .base_diffusion_limited import BaseModel
+
+
+class BaseFullDiffusionLimited(BaseModel):
+    """
+    Full submodel for diffusion-limited kinetics
+
+    Parameters
+    ----------
+    param :
+        model parameters
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
+
+
+    **Extends:** :class:`pybamm.interface.diffusion_limited.BaseModel`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def _get_diffusion_limited_current_density(self, variables):
+        param = self.param
+        if self.domain == "Negative":
+            N_ox_s_p = variables["Oxygen flux"].orphans[1]
+            N_ox_neg_sep_interface = N_ox_s_p[0]
+
+            j = -N_ox_neg_sep_interface / param.C_e / param.s_ox_Ox / param.l_n
+
+        return j
diff --git a/pybamm/models/submodels/interface/diffusion_limited/leading_diffusion_limited.py b/pybamm/models/submodels/interface/diffusion_limited/leading_diffusion_limited.py
new file mode 100644
index 000000000..f93d51d7e
--- /dev/null
+++ b/pybamm/models/submodels/interface/diffusion_limited/leading_diffusion_limited.py
@@ -0,0 +1,34 @@
+#
+# Leading-order diffusion limited kinetics
+#
+from .base_diffusion_limited import BaseModel
+
+
+class BaseLeadingOrderDiffusionLimited(BaseModel):
+    """
+    Leading-order submodel for diffusion-limited kinetics
+
+    Parameters
+    ----------
+    param :
+        model parameters
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
+
+
+    **Extends:** :class:`pybamm.interface.diffusion_limited.BaseModel`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def _get_diffusion_limited_current_density(self, variables):
+        if self.domain == "Negative":
+            j_p = variables[
+                "Positive electrode"
+                + self.reaction_name
+                + " interfacial current density"
+            ].orphans[0]
+            j = -self.param.l_p * j_p / self.param.l_n
+
+        return j
diff --git a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
index c027e9b38..dff0960c0 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -55,8 +55,12 @@ class BaseInverseButlerVolmer(BaseInterface):
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
         if (
-            "Negative electrode interfacial current density" in variables
-            and "Positive electrode interfacial current density" in variables
+            "Negative electrode" + self.reaction_name + " interfacial current density"
+            in variables
+            and "Positive electrode"
+            + self.reaction_name
+            + " interfacial current density"
+            in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/kinetics/__init__.py b/pybamm/models/submodels/interface/kinetics/__init__.py
index aacb13ef4..b3d47f69d 100644
--- a/pybamm/models/submodels/interface/kinetics/__init__.py
+++ b/pybamm/models/submodels/interface/kinetics/__init__.py
@@ -1,3 +1,4 @@
 from .base_kinetics import BaseModel
 from .butler_volmer import BaseButlerVolmer
 from .tafel import BaseForwardTafel, BaseBackwardTafel
+from .no_reaction import BaseNoReaction
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index 671cd42e1..bf54f1f9f 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -1,5 +1,5 @@
 #
-# Bulter volmer class
+# Base kinetics class
 #
 
 import pybamm
@@ -24,16 +24,6 @@ class BaseModel(BaseInterface):
     def __init__(self, param, domain):
         super().__init__(param, domain)
 
-    def _get_delta_phi_s(self, variables):
-        "Calculate delta_phi_s, and derived variables, using phi_s and phi_e"
-        phi_s = variables[self.domain + " electrode potential"]
-        phi_e = variables[self.domain + " electrolyte potential"]
-        delta_phi_s = phi_s - phi_e
-        variables.update(
-            self._get_standard_surface_potential_difference_variables(delta_phi_s)
-        )
-        return variables
-
     def get_coupled_variables(self, variables):
         # Calculate delta_phi_s from phi_s and phi_e if it isn't already known
         if self.domain + " electrode surface potential difference" not in variables:
@@ -48,11 +38,8 @@ class BaseModel(BaseInterface):
         # Get open-circuit potential variables and reaction overpotential
         ocp, dUdT = self._get_open_circuit_potential(variables)
         eta_r = delta_phi_s - ocp
-
-        if self.domain == "Negative":
-            ne = self.param.ne_n
-        elif self.domain == "Positive":
-            ne = self.param.ne_p
+        # Get number of electrons in reaction
+        ne = self._get_number_of_electrons_in_reaction()
 
         j = self._get_kinetics(j0, ne, eta_r)
         j_tot_av = self._get_average_total_interfacial_current_density(variables)
@@ -67,8 +54,12 @@ class BaseModel(BaseInterface):
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
         if (
-            "Negative electrode interfacial current density" in variables
-            and "Positive electrode interfacial current density" in variables
+            "Negative electrode" + self.reaction_name + " interfacial current density"
+            in variables
+            and "Positive electrode"
+            + self.reaction_name
+            + " interfacial current density"
+            in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/kinetics/no_reaction.py b/pybamm/models/submodels/interface/kinetics/no_reaction.py
new file mode 100644
index 000000000..c13f51321
--- /dev/null
+++ b/pybamm/models/submodels/interface/kinetics/no_reaction.py
@@ -0,0 +1,28 @@
+#
+# Bulter volmer class
+#
+
+import pybamm
+from .base_kinetics import BaseModel
+
+
+class BaseNoReaction(BaseModel):
+    """
+    Base submodel for when no reaction occurs
+
+    Parameters
+    ----------
+    param :
+        model parameters
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
+
+
+    **Extends:** :class:`pybamm.interface.kinetics.BaseKinetics`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def _get_kinetics(self, j0, ne, eta_r):
+        return pybamm.Scalar(0)
diff --git a/pybamm/models/submodels/interface/lead_acid.py b/pybamm/models/submodels/interface/lead_acid.py
index a540e2961..f4bb12abf 100644
--- a/pybamm/models/submodels/interface/lead_acid.py
+++ b/pybamm/models/submodels/interface/lead_acid.py
@@ -23,6 +23,7 @@ class BaseInterfaceLeadAcid(BaseInterface):
 
     def __init__(self, param, domain):
         super().__init__(param, domain)
+        self.reaction_name = ""  # empty reaction name, assumed to be the main reaction
 
     def _get_exchange_current_density(self, variables):
         """
@@ -55,19 +56,20 @@ class BaseInterfaceLeadAcid(BaseInterface):
 
     def _get_open_circuit_potential(self, variables):
         """
-        A private function to obtain the open circuit potential and
-        related standard variables.
+        A private function to obtain the open circuit potential and entropic change
 
         Parameters
         ----------
-        c_e : :class:`pybamm.Symbol`
-            The concentration in the electrolyte.
+        variables: dict
+            The variables in the full model.
 
         Returns
         -------
-        variables : dict
-            The variables dictionary including the open circuit potentials
-            and related standard variables.
+        ocp : :class:`pybamm.Symbol`
+            The open-circuit potential
+        dUdT : :class:`pybamm.Symbol`
+            The entropic change in open-circuit potential due to temperature
+
         """
 
         c_e = variables[self.domain + " electrolyte concentration"]
@@ -84,10 +86,17 @@ class BaseInterfaceLeadAcid(BaseInterface):
 
         return ocp, dUdT
 
+    def _get_number_of_electrons_in_reaction(self):
+        if self.domain == "Negative":
+            ne = self.param.ne_n
+        elif self.domain == "Positive":
+            ne = self.param.ne_p
+        return ne
+
 
 class ButlerVolmer(BaseInterfaceLeadAcid, kinetics.BaseButlerVolmer):
     """
-    Extends :class:`BaseInterfaceLeadIon` (for exchange-current density, etc) and
+    Extends :class:`BaseInterfaceLeadAcid` (for exchange-current density, etc) and
     :class:`kinetics.BaseButlerVolmer` (for kinetics)
     """
 
diff --git a/pybamm/models/submodels/interface/lead_acid_oxygen.py b/pybamm/models/submodels/interface/lead_acid_oxygen.py
new file mode 100644
index 000000000..c2a802f22
--- /dev/null
+++ b/pybamm/models/submodels/interface/lead_acid_oxygen.py
@@ -0,0 +1,126 @@
+#
+# Interface classes for oxygen reaction in lead-acid batteries
+#
+import pybamm
+from .base_interface import BaseInterface
+from . import kinetics, diffusion_limited
+
+
+class BaseInterfaceOxygenLeadAcid(BaseInterface):
+    """
+    Base lead-acid interface class
+
+    Parameters
+    ----------
+    param :
+        model parameters
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
+
+
+    **Extends:** :class:`pybamm.interface.BaseInterface`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+        self.reaction_name = " oxygen"
+
+    def _get_exchange_current_density(self, variables):
+        """
+        A private function to obtain the exchange current density for a lead acid
+        deposition reaction.
+
+        Parameters
+        ----------
+        variables: dict
+            The variables in the full model.
+
+        Returns
+        -------
+        j0 : :class: `pybamm.Symbol`
+            The exchange current density.
+        """
+        c_e = variables[self.domain + " electrolyte concentration"]
+        # If c_e was broadcast, take only the orphan
+        if isinstance(c_e, pybamm.Broadcast):
+            c_e = c_e.orphans[0]
+
+        if self.domain == "Negative":
+            j0 = pybamm.Scalar(0)
+        elif self.domain == "Positive":
+            j0 = self.param.j0_p_Ox_ref * c_e  # ** self.param.exponent_e_Ox
+
+        return j0
+
+    def _get_open_circuit_potential(self, variables):
+        """
+        A private function to obtain the open circuit potential and entropic change
+
+        Parameters
+        ----------
+        variables: dict
+            The variables in the full model.
+
+        Returns
+        -------
+        ocp : :class:`pybamm.Symbol`
+            The open-circuit potential
+        dUdT : :class:`pybamm.Symbol`
+            The entropic change in open-circuit potential due to temperature
+
+        """
+        if self.domain == "Negative":
+            ocp = self.param.U_n_Ox
+        elif self.domain == "Positive":
+            ocp = self.param.U_p_Ox
+
+        dUdT = pybamm.Scalar(0)
+
+        return ocp, dUdT
+
+    def _get_number_of_electrons_in_reaction(self):
+        return self.param.ne_Ox
+
+
+class ForwardTafel(BaseInterfaceOxygenLeadAcid, kinetics.BaseForwardTafel):
+    """
+    Extends :class:`BaseInterfaceOxygenLeadAcid` (for exchange-current density, etc) and
+    :class:`kinetics.BaseForwardTafel` (for kinetics)
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+
+class LeadingOrderDiffusionLimited(
+    BaseInterfaceOxygenLeadAcid, diffusion_limited.BaseLeadingOrderDiffusionLimited
+):
+    """
+    Extends :class:`BaseInterfaceOxygenLeadAcid` (for exchange-current density, etc) and
+    :class:`kinetics.BaseLeadingOrderDiffusionLimited` (for kinetics)
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+
+class FullDiffusionLimited(
+    BaseInterfaceOxygenLeadAcid, diffusion_limited.BaseFullDiffusionLimited
+):
+    """
+    Extends :class:`BaseInterfaceOxygenLeadAcid` (for exchange-current density, etc) and
+    :class:`kinetics.BaseFullDiffusionLimited` (for kinetics)
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+
+class NoReaction(BaseInterfaceOxygenLeadAcid, kinetics.BaseNoReaction):
+    """
+    Extends :class:`BaseInterfaceOxygenLeadAcid` (for exchange-current density, etc) and
+    :class:`kinetics.BaseNoReaction` (for kinetics)
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
diff --git a/pybamm/models/submodels/interface/lithium_ion.py b/pybamm/models/submodels/interface/lithium_ion.py
index 19277e48f..f17d1fd9b 100644
--- a/pybamm/models/submodels/interface/lithium_ion.py
+++ b/pybamm/models/submodels/interface/lithium_ion.py
@@ -22,6 +22,7 @@ class BaseInterfaceLithiumIon(BaseInterface):
 
     def __init__(self, param, domain):
         super().__init__(param, domain)
+        self.reaction_name = ""  # empty reaction name, assumed to be the main reaction
 
     def _get_exchange_current_density(self, variables):
         """
@@ -54,6 +55,22 @@ class BaseInterfaceLithiumIon(BaseInterface):
         return j0
 
     def _get_open_circuit_potential(self, variables):
+        """
+        A private function to obtain the open circuit potential and entropic change
+
+        Parameters
+        ----------
+        variables: dict
+            The variables in the full model.
+
+        Returns
+        -------
+        ocp : :class:`pybamm.Symbol`
+            The open-circuit potential
+        dUdT : :class:`pybamm.Symbol`
+            The entropic change in open-circuit potential due to temperature
+
+        """
         c_s_surf = variables[self.domain + " particle surface concentration"]
 
         if self.domain == "Negative":
@@ -66,6 +83,13 @@ class BaseInterfaceLithiumIon(BaseInterface):
 
         return ocp, dUdT
 
+    def _get_number_of_electrons_in_reaction(self):
+        if self.domain == "Negative":
+            ne = self.param.ne_n
+        elif self.domain == "Positive":
+            ne = self.param.ne_p
+        return ne
+
 
 class ButlerVolmer(BaseInterfaceLithiumIon, kinetics.BaseButlerVolmer):
     """
diff --git a/pybamm/models/submodels/oxygen_diffusion/__init__.py b/pybamm/models/submodels/oxygen_diffusion/__init__.py
new file mode 100644
index 000000000..a5161f93f
--- /dev/null
+++ b/pybamm/models/submodels/oxygen_diffusion/__init__.py
@@ -0,0 +1,4 @@
+from .base_oxygen_diffusion import BaseModel
+from .leading_oxygen_diffusion import LeadingOrder
+from .full_oxygen_diffusion import Full
+from .no_oxygen import NoOxygen
diff --git a/pybamm/models/submodels/oxygen_diffusion/base_oxygen_diffusion.py b/pybamm/models/submodels/oxygen_diffusion/base_oxygen_diffusion.py
new file mode 100644
index 000000000..d16ba8e5e
--- /dev/null
+++ b/pybamm/models/submodels/oxygen_diffusion/base_oxygen_diffusion.py
@@ -0,0 +1,84 @@
+#
+# Base class for oxygen diffusion
+#
+import pybamm
+
+
+class BaseModel(pybamm.BaseSubModel):
+    """Base class for conservation of mass of oxygen.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+
+
+    **Extends:** :class:`pybamm.BaseSubModel`
+    """
+
+    def __init__(self, param, reactions=None):
+        super().__init__(param, reactions=reactions)
+
+    def _get_standard_concentration_variables(self, c_ox):
+        """
+        A private function to obtain the standard variables which
+        can be derived from the concentration of oxygen.
+
+        Parameters
+        ----------
+        c_ox : :class:`pybamm.Symbol`
+            The concentration of oxygen.
+        c_ox_av : :class:`pybamm.Symbol`
+            The cell-averaged concentration of oxygen.
+
+        Returns
+        -------
+        variables : dict
+            The variables which can be derived from the concentration in the
+            oxygen.
+        """
+
+        c_ox_av = pybamm.average(c_ox)
+        c_ox_typ = self.param.c_ox_typ
+        c_ox_n, c_ox_s, c_ox_p = c_ox.orphans
+
+        variables = {
+            "Oxygen concentration": c_ox,
+            "Oxygen concentration [mol.m-3]": c_ox_typ * c_ox,
+            "Average oxygen concentration": c_ox_av,
+            "Average oxygen concentration [mol.m-3]": c_ox_typ * c_ox_av,
+            "Negative oxygen concentration": c_ox_n,
+            "Negative oxygen concentration [mol.m-3]": c_ox_typ * c_ox_n,
+            "Separator oxygen concentration": c_ox_s,
+            "Separator oxygen concentration [mol.m-3]": c_ox_typ * c_ox_s,
+            "Positive oxygen concentration": c_ox_p,
+            "Positive oxygen concentration [mol.m-3]": c_ox_typ * c_ox_p,
+        }
+
+        return variables
+
+    def _get_standard_flux_variables(self, N_ox):
+        """
+        A private function to obtain the standard variables which
+        can be derived from the mass flux of oxygen.
+
+        Parameters
+        ----------
+        N_ox : :class:`pybamm.Symbol`
+            The flux of oxygen.
+
+        Returns
+        -------
+        variables : dict
+            The variables which can be derived from the flux of oxygen.
+        """
+
+        param = self.param
+        flux_scale = param.curlyD_ox * param.c_ox_typ / param.L_x
+
+        variables = {
+            "Oxygen flux": N_ox,
+            "Oxygen flux [mol.m-2.s-1]": N_ox * flux_scale,
+        }
+
+        return variables
diff --git a/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py b/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py
new file mode 100644
index 000000000..0f60d32ab
--- /dev/null
+++ b/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py
@@ -0,0 +1,117 @@
+#
+# Class for oxygen diffusion
+#
+import pybamm
+
+from .base_oxygen_diffusion import BaseModel
+
+
+def separator_and_positive_only(variable):
+    """Return only the separator and positive electrode children
+
+    Parameters
+    ----------
+    variable : :class:`pybamm.Concatenation`
+        Concatenation of variables in negative, separator, positive
+
+    Returns
+    -------
+    :class:`pybamm.Concatenation`
+        Concatenation of variables in separator and positive only
+    """
+    _, var_s, var_p = variable.orphans
+    return pybamm.Concatenation(var_s, var_p)
+
+
+class Full(BaseModel):
+    """Class for conservation of mass of oxygen. (Full refers to unreduced by
+    asymptotic methods)
+    In this model, extremely fast oxygen kinetics in the negative electrode imposes
+    zero oxygen concentration there, and so the oxygen variable only lives in the
+    separator and positive electrode. The boundary condition at the negative electrode/
+    separator interface is homogeneous Dirichlet.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+
+
+    **Extends:** :class:`pybamm.oxygen.stefan_maxwell.diffusion.BaseModel`
+    """
+
+    def __init__(self, param, reactions):
+        super().__init__(param, reactions)
+
+    def get_fundamental_variables(self):
+        # Oxygen concentration (oxygen concentration is zero in the negative electrode)
+        c_ox_n = pybamm.Broadcast(0, "negative electrode")
+        c_ox_s = pybamm.Variable("Separator oxygen concentration", ["separator"])
+        c_ox_p = pybamm.Variable(
+            "Positive oxygen concentration", ["positive electrode"]
+        )
+        c_ox_s_p = pybamm.Concatenation(c_ox_s, c_ox_p)
+        variables = {"Separator and positive electrode oxygen concentration": c_ox_s_p}
+
+        c_ox = pybamm.Concatenation(c_ox_n, c_ox_s, c_ox_p)
+        variables.update(self._get_standard_concentration_variables(c_ox))
+
+        return variables
+
+    def get_coupled_variables(self, variables):
+
+        eps = separator_and_positive_only(variables["Porosity"])
+        c_ox = variables["Separator and positive electrode oxygen concentration"]
+        # TODO: allow charge and convection?
+        v_box = pybamm.Scalar(0)
+
+        param = self.param
+
+        N_ox_diffusion = -(eps ** param.b) * param.curlyD_ox * pybamm.grad(c_ox)
+
+        N_ox = N_ox_diffusion + c_ox * v_box
+        # Flux in the negative electrode is zero
+        N_ox = pybamm.Concatenation(pybamm.Broadcast(0, "negative electrode"), N_ox)
+
+        variables.update(self._get_standard_flux_variables(N_ox))
+
+        return variables
+
+    def set_rhs(self, variables):
+
+        param = self.param
+
+        eps = separator_and_positive_only(variables["Porosity"])
+        deps_dt = separator_and_positive_only(variables["Porosity change"])
+        c_ox = variables["Separator and positive electrode oxygen concentration"]
+        N_ox = variables["Oxygen flux"].orphans[1]
+
+        source_terms = sum(
+            pybamm.Concatenation(
+                pybamm.Broadcast(0, "separator"),
+                reaction["Positive"]["s_ox"] * variables[reaction["Positive"]["aj"]],
+            )
+            for reaction in self.reactions.values()
+        )
+
+        self.rhs = {
+            c_ox: (1 / eps)
+            * (-pybamm.div(N_ox) / param.C_e + source_terms - c_ox * deps_dt)
+        }
+
+    def set_boundary_conditions(self, variables):
+
+        c_ox = variables["Separator and positive electrode oxygen concentration"]
+
+        self.boundary_conditions = {
+            c_ox: {
+                "left": (pybamm.Scalar(0), "Dirichlet"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+
+    def set_initial_conditions(self, variables):
+
+        c_ox = variables["Separator and positive electrode oxygen concentration"]
+
+        self.initial_conditions = {c_ox: self.param.c_ox_init}
diff --git a/pybamm/models/submodels/oxygen_diffusion/leading_oxygen_diffusion.py b/pybamm/models/submodels/oxygen_diffusion/leading_oxygen_diffusion.py
new file mode 100644
index 000000000..0bd05e3fd
--- /dev/null
+++ b/pybamm/models/submodels/oxygen_diffusion/leading_oxygen_diffusion.py
@@ -0,0 +1,79 @@
+#
+# Class for leading-order oxygen diffusion
+#
+import pybamm
+
+from .base_oxygen_diffusion import BaseModel
+
+
+class LeadingOrder(BaseModel):
+    """Class for conservation of mass of oxygen. (Leading refers to leading order
+    of asymptotic reduction)
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    reactions : dict
+        Dictionary of reaction terms
+
+    **Extends:** :class:`pybamm.oxygen.stefan_maxwell.diffusion.BaseModel`
+    """
+
+    def __init__(self, param, reactions):
+        super().__init__(param, reactions)
+
+    def get_fundamental_variables(self):
+        c_ox_av = pybamm.Variable("Average oxygen concentration")
+        c_ox_n = pybamm.Broadcast(c_ox_av, ["negative electrode"])
+        c_ox_s = pybamm.Broadcast(c_ox_av, ["separator"])
+        c_ox_p = pybamm.Broadcast(c_ox_av, ["positive electrode"])
+        c_ox = pybamm.Concatenation(c_ox_n, c_ox_s, c_ox_p)
+
+        return self._get_standard_concentration_variables(c_ox)
+
+    def get_coupled_variables(self, variables):
+
+        N_ox = pybamm.Broadcast(
+            0, ["negative electrode", "separator", "positive electrode"]
+        )
+
+        variables.update(self._get_standard_flux_variables(N_ox))
+
+        return variables
+
+    def set_rhs(self, variables):
+
+        param = self.param
+
+        c_ox_av = variables["Average oxygen concentration"]
+
+        eps_n_av = variables["Average negative electrode porosity"]
+        eps_s_av = variables["Average separator porosity"]
+        eps_p_av = variables["Average positive electrode porosity"]
+
+        deps_n_dt_av = variables["Average negative electrode porosity change"]
+        deps_p_dt_av = variables["Average positive electrode porosity change"]
+
+        source_terms = sum(
+            param.l_n
+            * rxn["Negative"]["s_ox"]
+            * variables[rxn["Negative"]["aj"]].orphans[0]
+            + param.l_p
+            * rxn["Positive"]["s_ox"]
+            * variables[rxn["Positive"]["aj"]].orphans[0]
+            for rxn in self.reactions.values()
+        )
+
+        self.rhs = {
+            c_ox_av: 1
+            / (param.l_n * eps_n_av + param.l_s * eps_s_av + param.l_p * eps_p_av)
+            * (
+                source_terms
+                - c_ox_av * (param.l_n * deps_n_dt_av + param.l_p * deps_p_dt_av)
+            )
+        }
+
+    def set_initial_conditions(self, variables):
+        c_ox = variables["Average oxygen concentration"]
+        self.initial_conditions = {c_ox: self.param.c_ox_init}
diff --git a/pybamm/models/submodels/oxygen_diffusion/no_oxygen.py b/pybamm/models/submodels/oxygen_diffusion/no_oxygen.py
new file mode 100644
index 000000000..09eadf480
--- /dev/null
+++ b/pybamm/models/submodels/oxygen_diffusion/no_oxygen.py
@@ -0,0 +1,38 @@
+#
+# Class for when there is no oxygen
+#
+import pybamm
+
+from .base_oxygen_diffusion import BaseModel
+
+
+class NoOxygen(BaseModel):
+    """Class for when there is no oxygen
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+
+
+    **Extends:** :class:`pybamm.oxygen_diffusion.BaseModel`
+    """
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def get_fundamental_variables(self):
+        c_ox_n = pybamm.Broadcast(0, ["negative electrode"])
+        c_ox_s = pybamm.Broadcast(0, ["separator"])
+        c_ox_p = pybamm.Broadcast(0, ["positive electrode"])
+        c_ox = pybamm.Concatenation(c_ox_n, c_ox_s, c_ox_p)
+
+        variables = self._get_standard_concentration_variables(c_ox)
+
+        N_e = pybamm.Broadcast(
+            0, ["negative electrode", "separator", "positive electrode"]
+        )
+
+        variables.update(self._get_standard_flux_variables(N_e))
+
+        return variables
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 845723ecd..04c71d8da 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -60,11 +60,23 @@ class ParameterValues(dict):
         return {k: v for (k, v) in zip(df["Name [units]"], df["Value"])}
 
     def update(self, values):
+        # check parameter values
+        self.check_parameter_values(values)
+        # update
         for k, v in values.items():
             self[k] = v
         # reset processed symbols
         self._processed_symbols = {}
 
+    def check_parameter_values(self, values):
+        if "Typical current [A]" in values and values["Typical current [A]"] == 0:
+            raise ValueError(
+                """
+                "Typical current [A]" cannot be zero. A possible alternative is to set
+                "Current function" to `pybamm.GetConstantCurrent(current=0)` instead
+                """
+            )
+
     def process_model(self, model, processing="process"):
         """Assign parameter values to a model.
         Currently inplace, could be changed to return a new model.


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py
index 3903a8de5..18ea50844 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py
@@ -41,7 +41,7 @@ class TestLeadAcidNewmanTiedemann(unittest.TestCase):
         np.testing.assert_array_almost_equal(original, simp_and_python)
 
 
-class TestLeadAcidNewmanTiedemannCapacitance(unittest.TestCase):
+class TestLeadAcidNewmanTiedemannSurfaceForm(unittest.TestCase):
     def test_basic_processing_differential(self):
         options = {"surface form": "differential"}
         model = pybamm.lead_acid.NewmanTiedemann(options)
diff --git a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann_side_reactions.py b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann_side_reactions.py
new file mode 100644
index 000000000..b46f8626a
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann_side_reactions.py
@@ -0,0 +1,71 @@
+#
+# Tests for the lead-acid Newman-Tiedemann model
+#
+import pybamm
+import tests
+
+import unittest
+import numpy as np
+
+
+class TestLeadAcidNewmanTiedemannSideReactions(unittest.TestCase):
+    def test_basic_processing(self):
+        options = {"side reactions": ["oxygen"]}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all(skip_output_tests=True)
+
+    def test_basic_processing_differential(self):
+        options = {"side reactions": ["oxygen"], "surface form": "differential"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all(skip_output_tests=True)
+
+    @unittest.skipIf(pybamm.have_scikits_odes(), "scikits.odes not installed")
+    def test_basic_processing_algebraic(self):
+        options = {"side reactions": ["oxygen"], "surface form": "algebraic"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all(skip_output_tests=True)
+
+    def test_basic_processing_charge(self):
+        options = {"side reactions": ["oxygen"], "surface form": "differential"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        parameter_values = model.default_parameter_values
+        parameter_values.update(
+            {"Typical current [A]": -1, "Initial State of Charge": 0.5}
+        )
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)
+
+    def test_basic_processing_zero_current(self):
+        options = {"side reactions": ["oxygen"], "surface form": "differential"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        parameter_values = model.default_parameter_values
+        parameter_values.update(
+            {"Current function": pybamm.GetConstantCurrent(current=0)}
+        )
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)
+
+    def test_optimisations(self):
+        options = {"side reactions": ["oxygen"], "surface form": "differential"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        optimtest = tests.OptimisationsTest(model)
+
+        original = optimtest.evaluate_model()
+        simplified = optimtest.evaluate_model(simplify=True)
+        using_known_evals = optimtest.evaluate_model(use_known_evals=True)
+        simp_and_known = optimtest.evaluate_model(simplify=True, use_known_evals=True)
+        np.testing.assert_array_almost_equal(original, simplified, decimal=5)
+        np.testing.assert_array_almost_equal(original, using_known_evals)
+        np.testing.assert_array_almost_equal(original, simp_and_known, decimal=5)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_expression_tree/test_simplify.py b/tests/unit/test_expression_tree/test_simplify.py
index a02969d6c..8a9ba035c 100644
--- a/tests/unit/test_expression_tree/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_simplify.py
@@ -236,24 +236,24 @@ class TestSimplify(unittest.TestCase):
             sym.simplify()
 
         # A + A = 2A (#323)
-        a = pybamm.Variable("A")
+        a = pybamm.Parameter("A")
         expr = (a + a).simplify()
         self.assertIsInstance(expr, pybamm.Multiplication)
         self.assertIsInstance(expr.children[0], pybamm.Scalar)
         self.assertEqual(expr.children[0].evaluate(), 2)
-        self.assertIsInstance(expr.children[1], pybamm.Variable)
+        self.assertIsInstance(expr.children[1], pybamm.Parameter)
 
         expr = (a + a + a + a).simplify()
         self.assertIsInstance(expr, pybamm.Multiplication)
         self.assertIsInstance(expr.children[0], pybamm.Scalar)
         self.assertEqual(expr.children[0].evaluate(), 4)
-        self.assertIsInstance(expr.children[1], pybamm.Variable)
+        self.assertIsInstance(expr.children[1], pybamm.Parameter)
 
         expr = (a - a + a - a + a + a).simplify()
         self.assertIsInstance(expr, pybamm.Multiplication)
         self.assertIsInstance(expr.children[0], pybamm.Scalar)
         self.assertEqual(expr.children[0].evaluate(), 2)
-        self.assertIsInstance(expr.children[1], pybamm.Variable)
+        self.assertIsInstance(expr.children[1], pybamm.Parameter)
 
         # A - A = 0 (#323)
         expr = (a - a).simplify()
@@ -461,8 +461,9 @@ class TestSimplify(unittest.TestCase):
 
         # A - A = 0 (#323)
         expr = (a - a).simplify()
-        self.assertIsInstance(expr, pybamm.Scalar)
-        self.assertEqual(expr.evaluate(), 0)
+        self.assertIsInstance(expr, pybamm.Vector)
+        self.assertEqual(expr.shape, a.shape)
+        np.testing.assert_array_equal(expr.evaluate(), 0)
 
         # zero matrix
         m1 = pybamm.Matrix(np.zeros((300, 300)))
@@ -490,6 +491,20 @@ class TestSimplify(unittest.TestCase):
                 expr_simp.evaluate(y=np.ones(300)), -m2.evaluate()
             )
 
+    def test_matrix_divide_simplify(self):
+        m = pybamm.Matrix(np.random.rand(300, 200))
+        zero = pybamm.Scalar(0)
+
+        expr1 = (zero / m).simplify()
+        self.assertIsInstance(expr1, pybamm.Matrix)
+        self.assertEqual(expr1.shape, m.shape)
+        np.testing.assert_array_equal(expr1.evaluate(), 0)
+
+        expr2 = (m / zero).simplify()
+        self.assertIsInstance(expr2, pybamm.Matrix)
+        self.assertEqual(expr2.shape, m.shape)
+        np.testing.assert_array_equal(expr2.evaluate(), np.inf)
+
     def test_domain_concatenation_simplify(self):
         # create discretisation
         disc = get_discretisation_for_testing()
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index 4bc91d730..eee492be7 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -248,12 +248,10 @@ class TestSymbol(unittest.TestCase):
         param = pybamm.standard_parameters_lithium_ion
 
         one_n = pybamm.Broadcast(1, ["negative electrode"])
-        zero_s = pybamm.Broadcast(0, ["separator"])
         one_p = pybamm.Broadcast(1, ["positive electrode"])
 
-        j = pybamm.Concatenation(one_n, zero_s, one_p)
-
         zero_n = pybamm.Broadcast(0, ["negative electrode"])
+        zero_s = pybamm.Broadcast(0, ["separator"])
         zero_p = pybamm.Broadcast(0, ["positive electrode"])
 
         deps_dt = pybamm.Concatenation(zero_n, zero_s, zero_p)
@@ -263,18 +261,18 @@ class TestSymbol(unittest.TestCase):
         variables = {
             "Porosity": param.epsilon,
             "Porosity change": deps_dt,
-            "Interfacial current density": j,
             "Volume-averaged velocity": v_box,
+            "Negative electrode interfacial current density": one_n,
+            "Positive electrode interfacial current density": one_p,
         }
-
-        # c_e = pybamm.standard_variables.c_e
-        # variables = {"Electrolyte concentration": c_e}
-        # onen = pybamm.Broadcast(1, ["negative electrode"])
-        # onep = pybamm.Broadcast(1, ["positive electrode"])
-        # reactions = {
-        #     "main": {"neg": {"s": 1, "aj": onen}, "pos": {"s": 1, "aj": onep}}
-
-        model = pybamm.electrolyte.stefan_maxwell.diffusion.Full(param)
+        icd = " interfacial current density"
+        reactions = {
+            "main": {
+                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
+                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
+            }
+        }
+        model = pybamm.electrolyte.stefan_maxwell.diffusion.Full(param, reactions)
         variables.update(model.get_fundamental_variables())
         variables.update(model.get_coupled_variables(variables))
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
index be8e83742..46b2cebda 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
@@ -35,11 +35,16 @@ class TestLeadAcidLOQS(unittest.TestCase):
 
 
 class TestLeadAcidLOQSWithSideReactions(unittest.TestCase):
-    def test_well_posed(self):
+    def test_well_posed_differential(self):
         options = {"surface form": "differential", "side reactions": ["oxygen"]}
         model = pybamm.lead_acid.LOQS(options)
         model.check_well_posedness()
 
+    def test_well_posed_algebraic(self):
+        options = {"surface form": "algebraic", "side reactions": ["oxygen"]}
+        model = pybamm.lead_acid.LOQS(options)
+        model.check_well_posedness()
+
     def test_varying_surface_area(self):
         options = {
             "surface form": "differential",
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py
index d51c69721..1ff867da7 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_newman_tiedemann.py
@@ -21,7 +21,7 @@ class TestLeadAcidNewmanTiedemann(unittest.TestCase):
         self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
 
 
-class TestLeadAcidNewmanTiedemannCapacitance(unittest.TestCase):
+class TestLeadAcidNewmanTiedemannSurfaceForm(unittest.TestCase):
     def test_well_posed_differential(self):
         options = {"surface form": "differential"}
         model = pybamm.lead_acid.NewmanTiedemann(options)
@@ -42,6 +42,23 @@ class TestLeadAcidNewmanTiedemannCapacitance(unittest.TestCase):
         self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
 
 
+class TestLeadAcidNewmanTiedemannSideReactions(unittest.TestCase):
+    def test_well_posed(self):
+        options = {"side reactions": ["oxygen"]}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        model.check_well_posedness()
+
+    def test_well_posed_surface_form_differential(self):
+        options = {"side reactions": ["oxygen"], "surface form": "differential"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        model.check_well_posedness()
+
+    def test_well_posed_surface_form_algebraic(self):
+        options = {"side reactions": ["oxygen"], "surface form": "algebraic"}
+        model = pybamm.lead_acid.NewmanTiedemann(options)
+        model.check_well_posedness()
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py
index 2ec423107..d05645d45 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py
@@ -14,9 +14,21 @@ class TestFull(unittest.TestCase):
         variables = {
             "Porosity": a,
             "Electrolyte concentration": a,
-            "Interfacial current density": a,
+            "Negative electrode interfacial current density": pybamm.Broadcast(
+                a, "negative electrode"
+            ),
+            "Positive electrode interfacial current density": pybamm.Broadcast(
+                a, "positive electrode"
+            ),
         }
-        submodel = pybamm.electrolyte.stefan_maxwell.conductivity.Full(param)
+        icd = " interfacial current density"
+        reactions = {
+            "main": {
+                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
+                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
+            }
+        }
+        submodel = pybamm.electrolyte.stefan_maxwell.conductivity.Full(param, reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
index e3d8cc7f5..1d6f0b613 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
@@ -20,12 +20,19 @@ class TestFull(unittest.TestCase):
             "Negative electrolyte concentration": a_n,
             "Negative electrode interfacial current density": a_n,
         }
+        icd = " interfacial current density"
+        reactions = {
+            "main": {
+                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
+                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
+            }
+        }
 
         spf = pybamm.electrolyte.stefan_maxwell.conductivity.surface_potential_form
-        submodel = spf.FullAlgebraic(param, "Negative")
+        submodel = spf.FullAlgebraic(param, "Negative", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
-        submodel = spf.FullDifferential(param, "Negative")
+        submodel = spf.FullDifferential(param, "Negative", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
@@ -39,10 +46,10 @@ class TestFull(unittest.TestCase):
             "Positive electrolyte concentration": a_p,
             "Positive electrode interfacial current density": a_p,
         }
-        submodel = spf.FullAlgebraic(param, "Positive")
+        submodel = spf.FullAlgebraic(param, "Positive", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
-        submodel = spf.FullDifferential(param, "Positive")
+        submodel = spf.FullDifferential(param, "Positive", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py
index d6e85f0ae..378bf57a7 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py
@@ -21,12 +21,18 @@ class TestLeadingOrderModel(unittest.TestCase):
             "Negative electrode interfacial current density": a_n,
             "Average negative electrode total interfacial current density": a,
         }
-
+        icd = " interfacial current density"
+        reactions = {
+            "main": {
+                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
+                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
+            }
+        }
         spf = pybamm.electrolyte.stefan_maxwell.conductivity.surface_potential_form
-        submodel = spf.LeadingOrderAlgebraic(param, "Negative")
+        submodel = spf.LeadingOrderAlgebraic(param, "Negative", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
-        submodel = spf.LeadingOrderDifferential(param, "Negative")
+        submodel = spf.LeadingOrderDifferential(param, "Negative", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
@@ -41,10 +47,10 @@ class TestLeadingOrderModel(unittest.TestCase):
             "Positive electrode interfacial current density": a_p,
             "Average positive electrode total interfacial current density": a,
         }
-        submodel = spf.LeadingOrderAlgebraic(param, "Positive")
+        submodel = spf.LeadingOrderAlgebraic(param, "Positive", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
-        submodel = spf.LeadingOrderDifferential(param, "Positive")
+        submodel = spf.LeadingOrderDifferential(param, "Positive", reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py
index 4a6cb1313..e8b795583 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py
@@ -15,9 +15,22 @@ class TestFull(unittest.TestCase):
             "Porosity": a,
             "Porosity change": a,
             "Volume-averaged velocity": a,
-            "Interfacial current density": a,
+            "Electrolyte concentration": a,
+            "Negative electrode interfacial current density": pybamm.Broadcast(
+                a, "negative electrode"
+            ),
+            "Positive electrode interfacial current density": pybamm.Broadcast(
+                a, "positive electrode"
+            ),
         }
-        submodel = pybamm.electrolyte.stefan_maxwell.diffusion.Full(param)
+        icd = " interfacial current density"
+        reactions = {
+            "main": {
+                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
+                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
+            }
+        }
+        submodel = pybamm.electrolyte.stefan_maxwell.diffusion.Full(param, reactions)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_leading_stefan_maxwell_diffusion.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_leading_stefan_maxwell_diffusion.py
index e9959d035..f602e1078 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_leading_stefan_maxwell_diffusion.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_leading_stefan_maxwell_diffusion.py
@@ -12,13 +12,13 @@ class TestLeadingOrder(unittest.TestCase):
         param = pybamm.standard_parameters_lead_acid
         reactions = {
             "main": {
-                "neg": {
-                    "s_plus": param.s_n,
-                    "j": "Average negative electrode interfacial current density",
+                "Negative": {
+                    "s": param.s_n,
+                    "aj": "Negative electrode interfacial current density",
                 },
-                "pos": {
-                    "s_plus": param.s_p,
-                    "j": "Average positive electrode interfacial current density",
+                "Positive": {
+                    "s": param.s_p,
+                    "aj": "Positive electrode interfacial current density",
                 },
             }
         }
@@ -30,8 +30,12 @@ class TestLeadingOrder(unittest.TestCase):
             "Average negative electrode porosity change": a,
             "Average separator porosity change": a,
             "Average positive electrode porosity change": a,
-            "Average negative electrode interfacial current density": a,
-            "Average positive electrode interfacial current density": a,
+            "Negative electrode interfacial current density": pybamm.Broadcast(
+                a, "negative electrode"
+            ),
+            "Positive electrode interfacial current density": pybamm.Broadcast(
+                a, "positive electrode"
+            ),
         }
         submodel = pybamm.electrolyte.stefan_maxwell.diffusion.LeadingOrder(
             param, reactions
diff --git a/tests/unit/test_models/test_submodels/test_interface/test_diffusion_limited/test_base_diffusion_limited.py b/tests/unit/test_models/test_submodels/test_interface/test_diffusion_limited/test_base_diffusion_limited.py
new file mode 100644
index 000000000..d21c7d52e
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_interface/test_diffusion_limited/test_base_diffusion_limited.py
@@ -0,0 +1,27 @@
+#
+# Test base kinetics submodel
+#
+
+import pybamm
+import unittest
+
+
+class TestBaseModel(unittest.TestCase):
+    def test_not_implemented(self):
+        submodel = pybamm.interface.diffusion_limited.BaseModel(None, None)
+        with self.assertRaises(NotImplementedError):
+            submodel._get_exchange_current_density(None)
+        with self.assertRaises(NotImplementedError):
+            submodel._get_diffusion_limited_current_density(None)
+        with self.assertRaises(NotImplementedError):
+            submodel._get_open_circuit_potential(None)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/__init__.py b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_base_oxygen_diffusion.py b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_base_oxygen_diffusion.py
new file mode 100644
index 000000000..49adebab7
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_base_oxygen_diffusion.py
@@ -0,0 +1,24 @@
+#
+# Test base stefan maxwell diffusion submodel
+#
+
+import pybamm
+import tests
+import unittest
+
+
+class TestBaseModel(unittest.TestCase):
+    def test_public_functions(self):
+        variables = {"Oxygen concentration": pybamm.Scalar(0)}
+        submodel = pybamm.oxygen_diffusion.BaseModel(None)
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py
new file mode 100644
index 000000000..1086d49ac
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py
@@ -0,0 +1,56 @@
+#
+# Test full concentration submodel
+#
+
+import pybamm
+import tests
+import unittest
+
+
+class TestFull(unittest.TestCase):
+    def test_public_functions(self):
+        param = pybamm.standard_parameters_lead_acid
+        a = pybamm.Scalar(0)
+        variables = {
+            "Porosity": pybamm.Concatenation(
+                pybamm.Broadcast(a, "negative electrode"),
+                pybamm.Broadcast(a, "separator"),
+                pybamm.Broadcast(a, "positive electrode"),
+            ),
+            "Porosity change": pybamm.Concatenation(
+                pybamm.Broadcast(a, "negative electrode"),
+                pybamm.Broadcast(a, "separator"),
+                pybamm.Broadcast(a, "positive electrode"),
+            ),
+            "Volume-averaged velocity": a,
+            "Negative electrode interfacial current density": pybamm.Broadcast(
+                a, "negative electrode"
+            ),
+            "Positive electrode interfacial current density": pybamm.Broadcast(
+                a, "positive electrode"
+            ),
+        }
+        reactions = {
+            "main": {
+                "Negative": {
+                    "s_ox": param.s_ox_Ox,
+                    "aj": "Negative electrode interfacial current density",
+                },
+                "Positive": {
+                    "s_ox": param.s_ox_Ox,
+                    "aj": "Positive electrode interfacial current density",
+                },
+            }
+        }
+        submodel = pybamm.oxygen_diffusion.Full(param, reactions)
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_leading_oxygen_diffusion.py b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_leading_oxygen_diffusion.py
new file mode 100644
index 000000000..c44d4f02a
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_leading_oxygen_diffusion.py
@@ -0,0 +1,51 @@
+#
+# Test leading-order concentration submodel
+#
+
+import pybamm
+import tests
+import unittest
+
+
+class TestLeadingOrder(unittest.TestCase):
+    def test_public_functions(self):
+        param = pybamm.standard_parameters_lead_acid
+        reactions = {
+            "main": {
+                "Negative": {
+                    "s_ox": param.s_ox_Ox,
+                    "aj": "Negative electrode interfacial current density",
+                },
+                "Positive": {
+                    "s_ox": param.s_ox_Ox,
+                    "aj": "Positive electrode interfacial current density",
+                },
+            }
+        }
+        a = pybamm.Scalar(0)
+        variables = {
+            "Average negative electrode porosity": a,
+            "Average separator porosity": a,
+            "Average positive electrode porosity": a,
+            "Average negative electrode porosity change": a,
+            "Average separator porosity change": a,
+            "Average positive electrode porosity change": a,
+            "Negative electrode interfacial current density": pybamm.Broadcast(
+                a, "negative electrode"
+            ),
+            "Positive electrode interfacial current density": pybamm.Broadcast(
+                a, "positive electrode"
+            ),
+        }
+        submodel = pybamm.oxygen_diffusion.LeadingOrder(param, reactions)
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_no_oxygen.py b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_no_oxygen.py
new file mode 100644
index 000000000..9e52da2fd
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_no_oxygen.py
@@ -0,0 +1,24 @@
+#
+# Test constant concentration submodel
+#
+
+import pybamm
+import tests
+import unittest
+
+
+class TestConstantConcentration(unittest.TestCase):
+    def test_public_functions(self):
+        param = pybamm.standard_parameters_lead_acid
+        submodel = pybamm.oxygen_diffusion.NoOxygen(param)
+        std_tests = tests.StandardSubModelTests(submodel)
+        std_tests.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 91d768ffd..c656d72c3 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -41,6 +41,12 @@ class TestParameterValues(unittest.TestCase):
         self.assertEqual(param["Reference temperature [K]"], 294.85)
         self.assertEqual(param["Negative electrode width [m]"], 0.5)
 
+    def test_check_parameter_values(self):
+        # Can't provide a current density of 0, as this will cause a ZeroDivision error
+        bad_values = {"Typical current [A]": 0}
+        with self.assertRaises(ValueError):
+            pybamm.ParameterValues(bad_values)
+
     def test_process_symbol(self):
         parameter_values = pybamm.ParameterValues({"a": 1, "b": 2, "c": 3})
         # process parameter