From feb4d4f6cce98ba3d25cc225ffda2746386fc3cf Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Tue, 12 May 2020 16:42:06 -0400
Subject: [PATCH 01/15] #984 add sei resistance to inverse BV

---
 examples/scripts/compare_lithium_ion.py       | 10 +++--
 .../lithium-ion/seis/Example/parameters.csv   |  2 +-
 .../full_battery_models/base_battery_model.py |  1 +
 .../lithium_ion/base_lithium_ion_model.py     |  5 +++
 .../inverse_kinetics/inverse_butler_volmer.py | 10 +++--
 pybamm/models/submodels/sei/__init__.py       |  1 +
 pybamm/models/submodels/sei/constant_sei.py   | 38 +++++++++++++++++++
 pybamm/parameters/sei_parameters.py           |  4 +-
 8 files changed, 61 insertions(+), 10 deletions(-)
 create mode 100644 pybamm/models/submodels/sei/constant_sei.py

diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py
index 42cd480390..921fef776b 100644
--- a/examples/scripts/compare_lithium_ion.py
+++ b/examples/scripts/compare_lithium_ion.py
@@ -16,17 +16,19 @@
     pybamm.set_logging_level("INFO")
 
 # load models
-options = {"thermal": "lumped"}
+# options = {"thermal": "lumped"}
 models = [
-    pybamm.lithium_ion.SPM(options),
-    pybamm.lithium_ion.SPMe(options),
-    pybamm.lithium_ion.DFN(options),
+    pybamm.lithium_ion.SPM(),
+    pybamm.lithium_ion.SPM({"sei": "constant"}),
+    # pybamm.lithium_ion.SPMe(options),
+    # pybamm.lithium_ion.DFN(options),
 ]
 
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
 param["Current function [A]"] = 1
+param["SEI resistivity [Ohm.m]"] = 1e6
 
 for model in models:
     param.process_model(model)
diff --git a/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv b/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
index 8abc7ca2f4..638b8d8aa3 100644
--- a/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
+++ b/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
@@ -6,7 +6,7 @@ Inner SEI reaction proportion,0.5,,
 Inner SEI partial molar volume [m3.mol-1],3e-6, Guess,
 Outer SEI partial molar volume [m3.mol-1],3e-6, Guess,
 SEI reaction exchange current density [A.m-2],1.5E-7, Guess,
-SEI resistance per unit thickness [Ohm.m-1],1, Ramadass,
+SEI resistivity [Ohm.m],1, Ramadass,
 Outer SEI solvent diffusivity [m2.s-1],2.5E-22, Single paper,
 Bulk solvent concentration [mol.m-3],2.636E3, Ploehn paper,
 Ratio of inner and outer SEI exchange current densities,1, Assume same,
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 3ca269e8fc..b05be7d5dc 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -53,6 +53,7 @@ class BaseBatteryModel(pybamm.BaseModel):
                 Set the sei submodel to be used. Options are:
 
                 - None: :class:`pybamm.sei.NoSEI` (no SEI growth)
+                - "constant": :class:`pybamm.sei.Constant` (constant SEI thickness)
                 - "reaction limited": :class:`pybamm.sei.ReactionLimited`
                 - "solvent diffusion limited": \
                     :class:`pybamm.sei.SolventDiffusionLimited`
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 7d20682889..7d69befcc8 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
@@ -44,6 +44,11 @@ def set_sei_submodel(self):
                 self.param, "Negative electrode"
             )
 
+        if self.options["sei"] == "constant":
+            self.submodels["negative sei"] = pybamm.sei.ConstantSEI(
+                self.param, "Negative electrode"
+            )
+
         elif self.options["sei"] == "reaction limited":
             self.submodels["negative sei"] = pybamm.sei.ReactionLimited(
                 self.param, "Negative electrode"
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 4e305418d0..56cf1205cb 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -48,7 +48,9 @@ def get_coupled_variables(self, variables):
         else:
             T = variables[self.domain + " electrode temperature"]
 
-        eta_r = self._get_overpotential(j, j0, ne, T)
+        L_sei = variables["Total " + self.domain.lower() + " electrode sei thickness"]
+
+        eta_r = self._get_overpotential(j, j0, ne, T, L_sei)
         delta_phi = eta_r + ocp
 
         variables.update(self._get_standard_interfacial_current_variables(j))
@@ -79,5 +81,7 @@ def get_coupled_variables(self, variables):
 
         return variables
 
-    def _get_overpotential(self, j, j0, ne, T):
-        return (2 * (1 + self.param.Theta * T) / ne) * pybamm.arcsinh(j / (2 * j0))
+    def _get_overpotential(self, j, j0, ne, T, L_sei):
+        return (2 * (1 + self.param.Theta * T) / ne) * pybamm.arcsinh(
+            j / (2 * j0)
+        ) + j * L_sei * pybamm.sei_parameters.R_sei
diff --git a/pybamm/models/submodels/sei/__init__.py b/pybamm/models/submodels/sei/__init__.py
index b4162f63d5..c6c0c29db4 100644
--- a/pybamm/models/submodels/sei/__init__.py
+++ b/pybamm/models/submodels/sei/__init__.py
@@ -1,5 +1,6 @@
 from .base_sei import BaseModel
 from .no_sei import NoSEI
+from .constant_sei import ConstantSEI
 from .reaction_limited import ReactionLimited
 from .solvent_diffusion_limited import SolventDiffusionLimited
 from .electron_migration_limited import ElectronMigrationLimited
diff --git a/pybamm/models/submodels/sei/constant_sei.py b/pybamm/models/submodels/sei/constant_sei.py
new file mode 100644
index 0000000000..43b450bee0
--- /dev/null
+++ b/pybamm/models/submodels/sei/constant_sei.py
@@ -0,0 +1,38 @@
+#
+# Class for constant SEI thickness
+#
+import pybamm
+from .base_sei import BaseModel
+
+
+class ConstantSEI(BaseModel):
+    """Base class for no SEI.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    domain : str
+        The domain of the model either 'Negative' or 'Positive'
+
+    **Extends:** :class:`pybamm.particle.BaseParticle`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def get_fundamental_variables(self):
+        param = self.param
+
+        # Constant thicknesses
+        L_inner = pybamm.sei_parameters.L_inner_0
+        L_outer = pybamm.sei_parameters.L_outer_0
+        variables = self._get_standard_thickness_variables(L_inner, L_outer)
+
+        # Reactions
+        zero = pybamm.FullBroadcast(
+            pybamm.Scalar(0), self.domain.lower(), "current collector"
+        )
+        variables.update(self._get_standard_reaction_variables(zero, zero))
+
+        return variables
diff --git a/pybamm/parameters/sei_parameters.py b/pybamm/parameters/sei_parameters.py
index 6b2931b02e..ef1106f7b0 100644
--- a/pybamm/parameters/sei_parameters.py
+++ b/pybamm/parameters/sei_parameters.py
@@ -12,7 +12,7 @@
 
 m_sei_dimensional = pybamm.Parameter("SEI reaction exchange current density [A.m-2]")
 
-R_sei_dimensional = pybamm.Parameter("SEI resistance per unit thickness [Ohm.m-1]")
+R_sei_dimensional = pybamm.Parameter("SEI resistivity [Ohm.m]")
 
 D_sol_dimensional = pybamm.Parameter("Outer SEI solvent diffusivity [m2.s-1]")
 c_sol_dimensional = pybamm.Parameter("Bulk solvent concentration [mol.m-3]")
@@ -72,7 +72,7 @@
 
 U_inner_electron = F * U_inner_dimensional / R / T_ref
 
-R_sei = F * i_typ * R_sei_dimensional * L_sei_0_dim / (a_n * L_x)
+R_sei = F * i_typ * R_sei_dimensional * L_sei_0_dim / (a_n * L_x) / R / T_ref
 
 v_bar = V_bar_outer_dimensional / V_bar_inner_dimensional
 

From 4643732110f0509deaef7aa9bc4b5930d91ca2b3 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 13 May 2020 16:23:07 -0400
Subject: [PATCH 02/15] #984 reformat SEI to add film resistance and loss of
 lithium

---
 docs/source/models/submodels/index.rst        |   1 -
 .../models/submodels/interface/index.rst      |   1 +
 .../{ => interface}/sei/base_sei.rst          |   0
 .../sei/electron_migration_limited.rst        |   0
 .../submodels/{ => interface}/sei/index.rst   |   0
 .../sei/interstitial_diffusion_limited.rst    |   0
 .../submodels/{ => interface}/sei/no_sei.rst  |   0
 .../{ => interface}/sei/reaction_limited.rst  |   0
 .../sei/solvent_diffusion_limited.rst         |   0
 examples/scripts/compare_lithium_ion.py       |  10 +-
 examples/scripts/sei_growth.py                |  43 +++++---
 pybamm/__init__.py                            |   2 +-
 .../full_battery_models/base_battery_model.py |  93 +++++++++++++---
 .../lithium_ion/base_lithium_ion_model.py     |  18 ++-
 .../full_battery_models/lithium_ion/dfn.py    |   4 +-
 .../full_battery_models/lithium_ion/spm.py    |  16 ++-
 .../full_battery_models/lithium_ion/spme.py   |   4 +-
 .../submodels/interface/base_interface.py     |  12 +-
 .../interface/inverse_kinetics/__init__.py    |   2 +-
 .../inverse_kinetics/inverse_butler_volmer.py |  76 ++++++++++---
 .../interface/kinetics/base_kinetics.py       |  95 +++++++++++++++-
 .../interface/kinetics/butler_volmer.py       |   7 +-
 .../submodels/interface/kinetics/tafel.py     |   7 +-
 .../submodels/{ => interface}/sei/__init__.py |   0
 .../submodels/{ => interface}/sei/base_sei.py | 104 ++++++++----------
 .../submodels/interface/sei/constant_sei.py   |  52 +++++++++
 .../sei/electron_migration_limited.py         |  36 +++---
 .../sei/interstitial_diffusion_limited.py     |  38 ++++---
 .../submodels/{ => interface}/sei/no_sei.py   |  17 ++-
 .../interface/sei/reaction_limited.py         |  94 ++++++++++++++++
 .../sei/solvent_diffusion_limited.py          |  34 +++---
 pybamm/models/submodels/sei/constant_sei.py   |  38 -------
 .../models/submodels/sei/reaction_limited.py  |  75 -------------
 .../standard_parameters_lithium_ion.py        |   2 -
 .../test_base_battery_model.py                |  15 ++-
 .../test_base_lead_acid_model.py              |   4 +-
 36 files changed, 594 insertions(+), 306 deletions(-)
 rename docs/source/models/submodels/{ => interface}/sei/base_sei.rst (100%)
 rename docs/source/models/submodels/{ => interface}/sei/electron_migration_limited.rst (100%)
 rename docs/source/models/submodels/{ => interface}/sei/index.rst (100%)
 rename docs/source/models/submodels/{ => interface}/sei/interstitial_diffusion_limited.rst (100%)
 rename docs/source/models/submodels/{ => interface}/sei/no_sei.rst (100%)
 rename docs/source/models/submodels/{ => interface}/sei/reaction_limited.rst (100%)
 rename docs/source/models/submodels/{ => interface}/sei/solvent_diffusion_limited.rst (100%)
 rename pybamm/models/submodels/{ => interface}/sei/__init__.py (100%)
 rename pybamm/models/submodels/{ => interface}/sei/base_sei.py (59%)
 create mode 100644 pybamm/models/submodels/interface/sei/constant_sei.py
 rename pybamm/models/submodels/{ => interface}/sei/electron_migration_limited.py (57%)
 rename pybamm/models/submodels/{ => interface}/sei/interstitial_diffusion_limited.py (55%)
 rename pybamm/models/submodels/{ => interface}/sei/no_sei.py (52%)
 create mode 100644 pybamm/models/submodels/interface/sei/reaction_limited.py
 rename pybamm/models/submodels/{ => interface}/sei/solvent_diffusion_limited.py (57%)
 delete mode 100644 pybamm/models/submodels/sei/constant_sei.py
 delete mode 100644 pybamm/models/submodels/sei/reaction_limited.py

diff --git a/docs/source/models/submodels/index.rst b/docs/source/models/submodels/index.rst
index 28281889c5..a72bfa243f 100644
--- a/docs/source/models/submodels/index.rst
+++ b/docs/source/models/submodels/index.rst
@@ -17,4 +17,3 @@ Submodels
     porosity/index
     thermal/index
     tortuosity/index
-    sei/index
diff --git a/docs/source/models/submodels/interface/index.rst b/docs/source/models/submodels/interface/index.rst
index 21585d2841..c2fb250218 100644
--- a/docs/source/models/submodels/interface/index.rst
+++ b/docs/source/models/submodels/interface/index.rst
@@ -9,3 +9,4 @@ Interface
   inverse_kinetics/index
   first_order_kinetics/index
   diffusion_limited
+  sei/index
diff --git a/docs/source/models/submodels/sei/base_sei.rst b/docs/source/models/submodels/interface/sei/base_sei.rst
similarity index 100%
rename from docs/source/models/submodels/sei/base_sei.rst
rename to docs/source/models/submodels/interface/sei/base_sei.rst
diff --git a/docs/source/models/submodels/sei/electron_migration_limited.rst b/docs/source/models/submodels/interface/sei/electron_migration_limited.rst
similarity index 100%
rename from docs/source/models/submodels/sei/electron_migration_limited.rst
rename to docs/source/models/submodels/interface/sei/electron_migration_limited.rst
diff --git a/docs/source/models/submodels/sei/index.rst b/docs/source/models/submodels/interface/sei/index.rst
similarity index 100%
rename from docs/source/models/submodels/sei/index.rst
rename to docs/source/models/submodels/interface/sei/index.rst
diff --git a/docs/source/models/submodels/sei/interstitial_diffusion_limited.rst b/docs/source/models/submodels/interface/sei/interstitial_diffusion_limited.rst
similarity index 100%
rename from docs/source/models/submodels/sei/interstitial_diffusion_limited.rst
rename to docs/source/models/submodels/interface/sei/interstitial_diffusion_limited.rst
diff --git a/docs/source/models/submodels/sei/no_sei.rst b/docs/source/models/submodels/interface/sei/no_sei.rst
similarity index 100%
rename from docs/source/models/submodels/sei/no_sei.rst
rename to docs/source/models/submodels/interface/sei/no_sei.rst
diff --git a/docs/source/models/submodels/sei/reaction_limited.rst b/docs/source/models/submodels/interface/sei/reaction_limited.rst
similarity index 100%
rename from docs/source/models/submodels/sei/reaction_limited.rst
rename to docs/source/models/submodels/interface/sei/reaction_limited.rst
diff --git a/docs/source/models/submodels/sei/solvent_diffusion_limited.rst b/docs/source/models/submodels/interface/sei/solvent_diffusion_limited.rst
similarity index 100%
rename from docs/source/models/submodels/sei/solvent_diffusion_limited.rst
rename to docs/source/models/submodels/interface/sei/solvent_diffusion_limited.rst
diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py
index 921fef776b..c92dfe2aa5 100644
--- a/examples/scripts/compare_lithium_ion.py
+++ b/examples/scripts/compare_lithium_ion.py
@@ -18,17 +18,19 @@
 # load models
 # options = {"thermal": "lumped"}
 models = [
-    pybamm.lithium_ion.SPM(),
-    pybamm.lithium_ion.SPM({"sei": "constant"}),
+    # pybamm.lithium_ion.SPMe(),
+    # pybamm.lithium_ion.SPMe({"sei": "constant"}),
     # pybamm.lithium_ion.SPMe(options),
-    # pybamm.lithium_ion.DFN(options),
+    pybamm.lithium_ion.DFN(),
+    pybamm.lithium_ion.DFN({"sei": "constant", "sei film resistance": "average"}),
+    pybamm.lithium_ion.DFN({"sei": "constant", "sei film resistance": "distributed"}),
 ]
 
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
 param["Current function [A]"] = 1
-param["SEI resistivity [Ohm.m]"] = 1e6
+param["SEI resistivity [Ohm.m]"] = 1
 
 for model in models:
     param.process_model(model)
diff --git a/examples/scripts/sei_growth.py b/examples/scripts/sei_growth.py
index 5fa549b115..c7c9c4f088 100644
--- a/examples/scripts/sei_growth.py
+++ b/examples/scripts/sei_growth.py
@@ -1,29 +1,37 @@
 import pybamm as pb
 import numpy as np
 
-pb.set_logging_level("INFO")
+pb.set_logging_level("DEBUG")
 
-options = {"sei": "reaction limited"}
-model = pb.lithium_ion.DFN(options)
+options = {
+    "sei": "reaction limited",
+    "sei film resistance": None,
+    # "surface form": "algebraic",
+}
+models = [pb.lithium_ion.SPM(options), pb.lithium_ion.DFN(options)]
 
-parameter_values = model.default_parameter_values
+sims = []
+for model in models:
+    parameter_values = model.default_parameter_values
 
-parameter_values["Current function [A]"] = 0
+    parameter_values["Current function [A]"] = 0
 
-sim = pb.Simulation(model, parameter_values=parameter_values)
+    sim = pb.Simulation(model, parameter_values=parameter_values)
 
-solver = pb.CasadiSolver(mode="fast")
+    solver = pb.CasadiSolver(mode="fast")
 
-years = 3
-days = years * 365
-hours = days * 24
-minutes = hours * 60
-seconds = minutes * 60
+    years = 3
+    days = years * 365
+    hours = days * 24
+    minutes = hours * 60
+    seconds = minutes * 60
 
-t_eval = np.linspace(0, seconds, 100)
+    t_eval = np.linspace(0, seconds, 100)
 
-sim.solve(t_eval=t_eval, solver=solver)
-sim.plot(
+    sim.solve(t_eval=t_eval, solver=solver)
+    sims.append(sim)
+pb.dynamic_plot(
+    sims,
     [
         "Terminal voltage [V]",
         "Negative particle surface concentration",
@@ -42,5 +50,8 @@
             "X-averaged negative electrode sei interfacial current density [A.m-2]",
             "X-averaged negative electrode interfacial current density [A.m-2]",
         ],
-    ]
+        "Sum of x-averaged negative electrode interfacial current densities",
+        "Sum of negative electrode interfacial current densities",
+        "X-averaged electrolyte concentration",
+    ],
 )
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 44033ec2df..c8459eaf33 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -134,8 +134,8 @@ def version(formatted=False):
     porosity,
     thermal,
     tortuosity,
-    sei,
 )
+from .models.submodels.interface import sei
 
 #
 # Geometry
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index b05be7d5dc..cba817bba0 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -55,12 +55,38 @@ class BaseBatteryModel(pybamm.BaseModel):
                 - None: :class:`pybamm.sei.NoSEI` (no SEI growth)
                 - "constant": :class:`pybamm.sei.Constant` (constant SEI thickness)
                 - "reaction limited": :class:`pybamm.sei.ReactionLimited`
-                - "solvent diffusion limited": \
+                - "solvent-diffusion limited": \
                     :class:`pybamm.sei.SolventDiffusionLimited`
-                - "electron migration limited": \
+                - "electron-migration limited": \
                     :class:`pybamm.sei.ElectronMigrationLimited`
-                - "lithium interstitial diffusion limited": \
+                - "interstitial-diffusion limited": \
                     :class:`pybamm.sei.InterstitialDiffusionLimited`
+            * "sei film resistance" : str
+                Set the submodel for additional term in the overpotential due to SEI.
+                The default value is "None" if the "sei" option is "None", and
+                "distributed" otherwise. This is because the "distributed" model is more
+                complex than the model with no additional resistance, which adds
+                unnecessary complexity if there is no SEI in the first place
+
+                - None: no additional resistance\
+
+                    .. math::
+                        \\eta_r = \\frac{F}{RT} * (\\phi_s - \\phi_e - U)
+
+                - "distributed": properly included additional resistance term\
+
+                    .. math::
+                        \\eta_r = \\frac{F}{RT}
+                        * (\\phi_s - \\phi_e - U - R_{sei} * L_{sei} * j)
+
+                - "average": constant additional resistance term (approximation to the \
+                    true model). This model can give similar results to the \
+                    "distributed" case without needing to make j an algebraic state\
+
+                    .. math::
+                        \\eta_r = \\frac{F}{RT}
+                        * (\\phi_s - \\phi_e - U - R_{sei} * L_{sei} * \\frac{I}{aL})
+
 
     **Extends:** :class:`pybamm.BaseModel`
     """
@@ -158,18 +184,28 @@ def options(self, extra_options):
             "external submodels": [],
             "sei": None,
         }
+        # Change the default for SEI film resistance based on which sei option is
+        # provided
+        extra_options = extra_options or {}
+        sei_option = extra_options.get("sei", None)  # return None if option not given
+        if sei_option is None:
+            default_options["sei film resistance"] = None
+        else:
+            default_options["sei film resistance"] = "distributed"
+        # The "sei film resistance" option will still be overridden by extra_options if
+        # provided
+
         options = pybamm.FuzzyDict(default_options)
         # any extra options overwrite the default options
-        if extra_options is not None:
-            for name, opt in extra_options.items():
-                if name in default_options:
-                    options[name] = opt
-                else:
-                    raise pybamm.OptionError(
-                        "Option '{}' not recognised. Best matches are {}".format(
-                            name, options.get_best_matches(name)
-                        )
+        for name, opt in extra_options.items():
+            if name in default_options:
+                options[name] = opt
+            else:
+                raise pybamm.OptionError(
+                    "Option '{}' not recognised. Best matches are {}".format(
+                        name, options.get_best_matches(name)
                     )
+                )
 
         # Options that are incompatible with models
         if isinstance(self, pybamm.lithium_ion.BaseModel):
@@ -183,6 +219,8 @@ def options(self, extra_options):
                     "Lead-acid models can only have thermal "
                     "effects if dimensionality is 0."
                 )
+            if options["sei"] is not None or options["sei film resistance"] is not None:
+                raise pybamm.OptionError("Lead-acid models cannot have SEI formation")
 
         # Some standard checks to make sure options are compatible
         if not (
@@ -235,6 +273,26 @@ def options(self, extra_options):
             raise pybamm.OptionError(
                 "Unknown thermal model '{}'".format(options["thermal"])
             )
+        if options["sei"] not in [
+            None,
+            "constant",
+            "reaction limited",
+            "solvent-diffusion limited",
+            "electron-migration limited",
+            "interstitial-diffusion limited",
+        ]:
+            raise pybamm.OptionError("Unknown sei model '{}'".format(options["sei"]))
+        if options["sei film resistance"] not in [
+            None,
+            "distributed",
+            "average",
+        ]:
+            raise pybamm.OptionError(
+                "Unknown sei film resistance model '{}'".format(
+                    options["sei film resistance"]
+                )
+            )
+
         if options["dimensionality"] == 0:
             if options["current collector"] not in [
                 "uniform",
@@ -253,8 +311,8 @@ def options(self, extra_options):
 
         if options["thermal"] == "x-lumped" and options["dimensionality"] == 1:
             warnings.warn(
-                "1+1D Thermal models are only valid if both tabs are"
-                + "placed at the top of the cell."
+                "1+1D Thermal models are only valid if both tabs are "
+                "placed at the top of the cell."
             )
 
         self._options = options
@@ -295,6 +353,9 @@ def set_standard_output_variables(self):
             )
 
         # Initialize "total reaction" variables
+        # These will get populated by the "get_coupled_variables" methods, and then used
+        # later by "set_rhs" or "set_algebraic", which ensures that we always have
+        # added all the necessary variables by the time the sum is used
         self.variables.update(
             {
                 "Sum of electrolyte reaction source terms": 0,
@@ -368,10 +429,10 @@ def build_coupled_variables(self):
                         if len(submodels) == 1 or count == 100:
                             # no more submodels to try
                             raise pybamm.ModelError(
-                                "Submodel '{}' requires the variable {}, ".format(
+                                "Missing variable for submodel '{}': {}.\n".format(
                                     submodel_name, key
                                 )
-                                + "but it cannot be found. Check the selected "
+                                + "Check the selected "
                                 "submodels provide all of the required variables."
                             )
                         else:
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 7d69befcc8..5176da2174 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
@@ -40,39 +40,35 @@ def set_sei_submodel(self):
 
         # negative electrode SEI
         if self.options["sei"] is None:
-            self.submodels["negative sei"] = pybamm.sei.NoSEI(
-                self.param, "Negative electrode"
-            )
+            self.submodels["negative sei"] = pybamm.sei.NoSEI(self.param, "Negative")
 
         if self.options["sei"] == "constant":
             self.submodels["negative sei"] = pybamm.sei.ConstantSEI(
-                self.param, "Negative electrode"
+                self.param, "Negative"
             )
 
         elif self.options["sei"] == "reaction limited":
             self.submodels["negative sei"] = pybamm.sei.ReactionLimited(
-                self.param, "Negative electrode"
+                self.param, "Negative"
             )
 
         elif self.options["sei"] == "solvent-diffusion limited":
             self.submodels["negative sei"] = pybamm.sei.SolventDiffusionLimited(
-                self.param, "Negative electrode"
+                self.param, "Negative"
             )
 
         elif self.options["sei"] == "electron-migration limited":
             self.submodels["negative sei"] = pybamm.sei.ElectronMigrationLimited(
-                self.param, "Negative electrode"
+                self.param, "Negative"
             )
 
         elif self.options["sei"] == "interstitial-diffusion limited":
             self.submodels["negative sei"] = pybamm.sei.InterstitialDiffusionLimited(
-                self.param, "Negative electrode"
+                self.param, "Negative"
             )
 
         # positive electrode
-        self.submodels["positive sei"] = pybamm.sei.NoSEI(
-            self.param, "Positive electrode"
-        )
+        self.submodels["positive sei"] = pybamm.sei.NoSEI(self.param, "Positive")
 
     def set_other_reaction_submodels_to_zero(self):
         self.submodels["negative oxygen interface"] = pybamm.interface.NoReaction(
diff --git a/pybamm/models/full_battery_models/lithium_ion/dfn.py b/pybamm/models/full_battery_models/lithium_ion/dfn.py
index 21245ddd20..85407406ed 100644
--- a/pybamm/models/full_battery_models/lithium_ion/dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/dfn.py
@@ -67,10 +67,10 @@ def set_convection_submodel(self):
     def set_interfacial_submodel(self):
 
         self.submodels["negative interface"] = pybamm.interface.ButlerVolmer(
-            self.param, "Negative", "lithium-ion main"
+            self.param, "Negative", "lithium-ion main", self.options
         )
         self.submodels["positive interface"] = pybamm.interface.ButlerVolmer(
-            self.param, "Positive", "lithium-ion main"
+            self.param, "Positive", "lithium-ion main", self.options
         )
 
     def set_particle_submodel(self):
diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py
index 8431f602c5..99178a32cf 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spm.py
@@ -68,18 +68,28 @@ def set_interfacial_submodel(self):
 
         if self.options["surface form"] is False:
             self.submodels["negative interface"] = pybamm.interface.InverseButlerVolmer(
-                self.param, "Negative", "lithium-ion main"
+                self.param, "Negative", "lithium-ion main", self.options
             )
             self.submodels["positive interface"] = pybamm.interface.InverseButlerVolmer(
+                self.param, "Positive", "lithium-ion main", self.options
+            )
+            self.submodels[
+                "negative interface current"
+            ] = pybamm.interface.CurrentForInverseButlerVolmer(
+                self.param, "Negative", "lithium-ion main"
+            )
+            self.submodels[
+                "positive interface current"
+            ] = pybamm.interface.CurrentForInverseButlerVolmer(
                 self.param, "Positive", "lithium-ion main"
             )
         else:
             self.submodels["negative interface"] = pybamm.interface.ButlerVolmer(
-                self.param, "Negative", "lithium-ion main"
+                self.param, "Negative", "lithium-ion main", self.options
             )
 
             self.submodels["positive interface"] = pybamm.interface.ButlerVolmer(
-                self.param, "Positive", "lithium-ion main"
+                self.param, "Positive", "lithium-ion main", self.options
             )
 
     def set_particle_submodel(self):
diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py
index 04e754e111..e033499988 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spme.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spme.py
@@ -78,10 +78,10 @@ def set_tortuosity_submodels(self):
     def set_interfacial_submodel(self):
 
         self.submodels["negative interface"] = pybamm.interface.InverseButlerVolmer(
-            self.param, "Negative", "lithium-ion main"
+            self.param, "Negative", "lithium-ion main", self.options
         )
         self.submodels["positive interface"] = pybamm.interface.InverseButlerVolmer(
-            self.param, "Positive", "lithium-ion main"
+            self.param, "Positive", "lithium-ion main", self.options
         )
 
     def set_particle_submodel(self):
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index 3f5aa35d9c..a6ed54dbd5 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -23,7 +23,6 @@ class BaseInterface(pybamm.BaseSubModel):
 
     def __init__(self, param, domain, reaction):
         super().__init__(param, domain)
-        self.reaction = reaction
         if reaction == "lithium-ion main":
             self.reaction_name = ""  # empty reaction name for the main reaction
         elif reaction == "lead-acid main":
@@ -32,6 +31,11 @@ def __init__(self, param, domain, reaction):
             self.reaction_name = " oxygen"
         elif reaction == "lithium-ion oxygen":
             self.reaction_name = " oxygen"
+        elif reaction == "sei":
+            self.reaction_name = " sei"
+        else:
+            raise ValueError("Reaction name '{}' not recognized".format(reaction))
+        self.reaction = reaction
 
     def _get_exchange_current_density(self, variables):
         """
@@ -168,7 +172,9 @@ def _get_number_of_electrons_in_reaction(self):
 
     def _get_electrolyte_reaction_signed_stoichiometry(self):
         "Returns the number of electrons in the reaction"
-        if self.reaction == "lithium-ion main":
+        if self.reaction in ["lithium-ion main", "sei"]:
+            # Both the main reaction current contribute to the electrolyte reaction
+            # current
             return pybamm.Scalar(1), pybamm.Scalar(1)
         elif self.reaction == "lead-acid main":
             return self.param.s_plus_n_S, self.param.s_plus_p_S
@@ -281,7 +287,7 @@ def _get_standard_total_interfacial_current_variables(self, j_tot_av):
 
     def _get_standard_whole_cell_interfacial_current_variables(self, variables):
         """
-        Get variables associated with interfacial current over theh whole cell domain
+        Get variables associated with interfacial current over the whole cell domain
         This function also automatically increments the "total source term" variables
         """
         i_typ = self.param.i_typ
diff --git a/pybamm/models/submodels/interface/inverse_kinetics/__init__.py b/pybamm/models/submodels/interface/inverse_kinetics/__init__.py
index 1283fde47b..3c81882350 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/__init__.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/__init__.py
@@ -1 +1 @@
-from .inverse_butler_volmer import InverseButlerVolmer
+from .inverse_butler_volmer import InverseButlerVolmer, CurrentForInverseButlerVolmer
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 56cf1205cb..da59aa9740 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -7,7 +7,7 @@
 
 class InverseButlerVolmer(BaseInterface):
     """
-    A base submodel that implements the inverted form of the Butler-Volmer relation to
+    A submodel that implements the inverted form of the Butler-Volmer relation to
     solve for the reaction overpotential.
 
     Parameters
@@ -19,13 +19,19 @@ class InverseButlerVolmer(BaseInterface):
         in which case j.domain is used.
     reaction : str
         The name of the reaction being implemented
+    options: dict
+        A dictionary of options to be passed to the model. In this case "sei film
+        resistance" is the important option. See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.interface.kinetics.ButlerVolmer`
 
     """
 
-    def __init__(self, param, domain, reaction):
+    def __init__(self, param, domain, reaction, options=None):
         super().__init__(param, domain, reaction)
+        if options is None:
+            options = {"sei film resistance": None}
+        self.options = options
 
     def get_coupled_variables(self, variables):
         ocp, dUdT = self._get_open_circuit_potential(variables)
@@ -34,9 +40,11 @@ def get_coupled_variables(self, variables):
         j_tot_av = self._get_average_total_interfacial_current_density(variables)
         # Broadcast to match j0's domain
         if j0.domain in [[], ["current collector"]]:
-            j = j_tot_av
+            j_tot = j_tot_av
         else:
-            j = pybamm.PrimaryBroadcast(j_tot_av, [self.domain.lower() + " electrode"])
+            j_tot = pybamm.PrimaryBroadcast(
+                j_tot_av, [self.domain.lower() + " electrode"]
+            )
 
         if self.domain == "Negative":
             ne = self.param.ne_n
@@ -48,12 +56,21 @@ def get_coupled_variables(self, variables):
         else:
             T = variables[self.domain + " electrode temperature"]
 
-        L_sei = variables["Total " + self.domain.lower() + " electrode sei thickness"]
+        # eta_r is the overpotential from inverting Butler-Volmer, regardless of any
+        # additional SEI resistance. What changes is how delta_phi is defined in terms
+        # of eta_r
+        eta_r = self._get_overpotential(j_tot, j0, ne, T)
 
-        eta_r = self._get_overpotential(j, j0, ne, T, L_sei)
-        delta_phi = eta_r + ocp
+        # With SEI resistance (distributed and averaged have the same effect here)
+        if self.options["sei film resistance"] is not None:
+            L_sei = variables[
+                "Total " + self.domain.lower() + " electrode sei thickness"
+            ]
+            delta_phi = eta_r + ocp + j_tot * L_sei * pybamm.sei_parameters.R_sei
+        # Without SEI resistance
+        else:
+            delta_phi = eta_r + ocp
 
-        variables.update(self._get_standard_interfacial_current_variables(j))
         variables.update(
             self._get_standard_total_interfacial_current_variables(j_tot_av)
         )
@@ -64,6 +81,44 @@ def get_coupled_variables(self, variables):
         )
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
+        return variables
+
+    def _get_overpotential(self, j, j0, ne, T):
+        return (2 * (1 + self.param.Theta * T) / ne) * pybamm.arcsinh(j / (2 * j0))
+
+
+class CurrentForInverseButlerVolmer(BaseInterface):
+    """
+    Submodel for the current associated with the inverse Butler-Volmer formulation. 
+
+    Parameters
+    ----------
+    param
+        Model parameters
+    domain : iter of str, optional
+        The domain(s) in which to compute the interfacial current. Default is None,
+        in which case j.domain is used.
+    reaction : str
+        The name of the reaction being implemented
+
+    **Extends:** :class:`pybamm.interface.kinetics.ButlerVolmer`
+
+    """
+
+    def __init__(self, param, domain, reaction):
+        super().__init__(param, domain, reaction)
+
+    def get_coupled_variables(self, variables):
+        j_tot = variables[
+            "X-averaged "
+            + self.domain.lower()
+            + " electrode total interfacial current density"
+        ]
+        j_sei = variables[self.domain + " electrode sei interfacial current density"]
+        j = j_tot - j_sei
+
+        variables.update(self._get_standard_interfacial_current_variables(j))
+
         if (
             "Negative electrode" + self.reaction_name + " interfacial current density"
             in variables
@@ -80,8 +135,3 @@ def get_coupled_variables(self, variables):
             )
 
         return variables
-
-    def _get_overpotential(self, j, j0, ne, T, L_sei):
-        return (2 * (1 + self.param.Theta * T) / ne) * pybamm.arcsinh(
-            j / (2 * j0)
-        ) + j * L_sei * pybamm.sei_parameters.R_sei
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index 6a32b646fe..1f67c28629 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -18,12 +18,34 @@ class BaseKinetics(BaseInterface):
         The domain to implement the model, either: 'Negative' or 'Positive'.
     reaction : str
         The name of the reaction being implemented
+    options: dict
+        A dictionary of options to be passed to the model. In this case "sei film
+        resistance" is the important option. See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.interface.BaseInterface`
     """
 
-    def __init__(self, param, domain, reaction):
+    def __init__(self, param, domain, reaction, options=None):
         super().__init__(param, domain, reaction)
+        if options is None:
+            options = {"sei film resistance": None}
+        self.options = options
+
+    def get_fundamental_variables(self):
+        if (
+            self.options["sei film resistance"] == "distributed"
+            and "main" in self.reaction
+        ):
+            j = pybamm.Variable(
+                self.domain + " electrode interfacial current density",
+                domain=self.domain.lower() + " electrode",
+                auxiliary_domains={"secondary": "current collector"},
+            )
+
+            variables = self._get_standard_interfacial_current_variables(j)
+            return variables
+        else:
+            return {}
 
     def get_coupled_variables(self, variables):
         # Calculate delta_phi from phi_s and phi_e if it isn't already known
@@ -39,6 +61,24 @@ def get_coupled_variables(self, variables):
         # Get open-circuit potential variables and reaction overpotential
         ocp, dUdT = self._get_open_circuit_potential(variables)
         eta_r = delta_phi - ocp
+
+        # Get average interfacial current density
+        j_tot_av = self._get_average_total_interfacial_current_density(variables)
+        # j = j_tot_av + (j - pybamm.x_average(j))  # enforce true average
+
+        # Add SEI resistance
+        if self.options["sei film resistance"] == "distributed":
+            L_sei = variables[
+                "Total " + self.domain.lower() + " electrode sei thickness"
+            ]
+            j = variables[self.domain + " electrode interfacial current density"]
+            eta_r -= j * L_sei * pybamm.sei_parameters.R_sei
+        elif self.options["sei film resistance"] == "average":
+            L_sei = variables[
+                "Total " + self.domain.lower() + " electrode sei thickness"
+            ]
+            eta_r -= j_tot_av * L_sei * pybamm.sei_parameters.R_sei
+
         # Get number of electrons in reaction
         ne = self._get_number_of_electrons_in_reaction()
         # Get kinetics. Note: T must have the same domain as j0 and eta_r
@@ -46,12 +86,21 @@ def get_coupled_variables(self, variables):
             T = variables["X-averaged cell temperature"]
         else:
             T = variables[self.domain + " electrode temperature"]
+
+        # Update j, except in the "distributed SEI resistance" model, where j will be
+        # found by solving an algebraic equation
+        # (In the "distributed SEI resistance" model, we have already defined j)
         j = self._get_kinetics(j0, ne, eta_r, T)
-        # Get average interfacial current density
-        j_tot_av = self._get_average_total_interfacial_current_density(variables)
-        # j = j_tot_av + (j - pybamm.x_average(j))  # enforce true average
+        if (
+            self.options["sei film resistance"] == "distributed"
+            and "main" in self.reaction
+        ):
+            variables.update(
+                {self.domain + " electrode" + self.reaction_name + " reaction term": j}
+            )
+        else:
+            variables.update(self._get_standard_interfacial_current_variables(j))
 
-        variables.update(self._get_standard_interfacial_current_variables(j))
         variables.update(
             self._get_standard_total_interfacial_current_variables(j_tot_av)
         )
@@ -76,6 +125,42 @@ def get_coupled_variables(self, variables):
 
         return variables
 
+    def set_algebraic(self, variables):
+        if (
+            self.options["sei film resistance"] == "distributed"
+            and "main" in self.reaction
+        ):
+            j = variables[
+                self.domain
+                + " electrode"
+                + self.reaction_name
+                + " interfacial current density"
+            ]
+            reaction = variables[
+                self.domain + " electrode" + self.reaction_name + " reaction term"
+            ]
+            # Algebraic equation to set j equal to the reaction term
+            self.algebraic[j] = j - reaction
+
+    def set_initial_conditions(self, variables):
+        if (
+            self.options["sei film resistance"] == "distributed"
+            and "main" in self.reaction
+        ):
+            param = self.param
+            j = variables[
+                self.domain
+                + " electrode"
+                + self.reaction_name
+                + " interfacial current density"
+            ]
+            if self.domain == "Negative":
+                j_av_init = param.current_with_time / param.l_n
+            elif self.domain == "Positive":
+                j_av_init = -param.current_with_time / param.l_p
+
+            self.initial_conditions[j] = j_av_init
+
     def _get_dj_dc(self, variables):
         """
         Default to calculate derivative of interfacial current density with respect to
diff --git a/pybamm/models/submodels/interface/kinetics/butler_volmer.py b/pybamm/models/submodels/interface/kinetics/butler_volmer.py
index 9a6480c78d..4acc05cac6 100644
--- a/pybamm/models/submodels/interface/kinetics/butler_volmer.py
+++ b/pybamm/models/submodels/interface/kinetics/butler_volmer.py
@@ -21,12 +21,15 @@ class ButlerVolmer(BaseKinetics):
         The domain to implement the model, either: 'Negative' or 'Positive'.
     reaction : str
         The name of the reaction being implemented
+    options: dict
+        A dictionary of options to be passed to the model. In this case "sei film
+        resistance" is the important option. See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.interface.kinetics.BaseKinetics`
     """
 
-    def __init__(self, param, domain, reaction):
-        super().__init__(param, domain, reaction)
+    def __init__(self, param, domain, reaction, options=None):
+        super().__init__(param, domain, reaction, options)
 
     def _get_kinetics(self, j0, ne, eta_r, T):
         prefactor = ne / (2 * (1 + self.param.Theta * T))
diff --git a/pybamm/models/submodels/interface/kinetics/tafel.py b/pybamm/models/submodels/interface/kinetics/tafel.py
index 4aff0c6a75..01e9c3f8ba 100644
--- a/pybamm/models/submodels/interface/kinetics/tafel.py
+++ b/pybamm/models/submodels/interface/kinetics/tafel.py
@@ -20,12 +20,15 @@ class ForwardTafel(BaseKinetics):
         The domain to implement the model, either: 'Negative' or 'Positive'.
     reaction : str
         The name of the reaction being implemented
+    options: dict
+        A dictionary of options to be passed to the model. In this case "sei film
+        resistance" is the important option. See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.interface.kinetics.BaseKinetics`
     """
 
-    def __init__(self, param, domain, reaction):
-        super().__init__(param, domain, reaction)
+    def __init__(self, param, domain, reaction, options=None):
+        super().__init__(param, domain, reaction, options)
 
     def _get_kinetics(self, j0, ne, eta_r, T):
         return j0 * pybamm.exp((ne / (2 * (1 + self.param.Theta * T))) * eta_r)
diff --git a/pybamm/models/submodels/sei/__init__.py b/pybamm/models/submodels/interface/sei/__init__.py
similarity index 100%
rename from pybamm/models/submodels/sei/__init__.py
rename to pybamm/models/submodels/interface/sei/__init__.py
diff --git a/pybamm/models/submodels/sei/base_sei.py b/pybamm/models/submodels/interface/sei/base_sei.py
similarity index 59%
rename from pybamm/models/submodels/sei/base_sei.py
rename to pybamm/models/submodels/interface/sei/base_sei.py
index 3299ae11f0..aee9408d06 100644
--- a/pybamm/models/submodels/sei/base_sei.py
+++ b/pybamm/models/submodels/interface/sei/base_sei.py
@@ -2,24 +2,25 @@
 # Base class for SEI models.
 #
 import pybamm
+from ..base_interface import BaseInterface
 
 
-class BaseModel(pybamm.BaseSubModel):
+class BaseModel(BaseInterface):
     """Base class for SEI models.
 
     Parameters
     ----------
     param : parameter class
         The parameters to use for this submodel
-    reactions : dict, optional
-        Dictionary of reaction terms
+    domain : str
+        The domain to implement the model, either: 'Negative' or 'Positive'.
 
-    **Extends:** :class:`pybamm.BaseSubModel`
+    **Extends:** :class:`pybamm.interface.BaseInterface`
     """
 
     def __init__(self, param, domain):
-        self.domain = domain
-        super().__init__(param)
+        reaction = "sei"
+        super().__init__(param, domain, reaction)
 
     def _get_standard_thickness_variables(self, L_inner, L_outer):
         """
@@ -69,44 +70,32 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
 
         Q_sei = n_SEI_av * self.param.L_n * self.param.L_y * self.param.L_z
 
+        domain = self.domain.lower() + " electrode"
+
         variables = {
-            "Inner " + self.domain.lower() + " sei thickness": L_inner,
-            "Inner " + self.domain.lower() + " sei thickness [m]": L_inner * L_scale,
-            "X-averaged inner " + self.domain.lower() + " sei thickness": L_inner_av,
-            "X-averaged inner "
-            + self.domain.lower()
-            + " sei thickness [m]": L_inner_av * L_scale,
-            "Outer " + self.domain.lower() + " sei thickness": L_outer,
-            "Outer " + self.domain.lower() + " sei thickness [m]": L_outer * L_scale,
-            "X-averaged outer " + self.domain.lower() + " sei thickness": L_outer_av,
-            "X-averaged outer "
-            + self.domain.lower()
-            + " sei thickness [m]": L_outer_av * L_scale,
-            "Total " + self.domain.lower() + " sei thickness": L_sei,
-            "Total " + self.domain.lower() + " sei thickness [m]": L_sei * L_scale,
-            "X-averaged total " + self.domain.lower() + " sei thickness": L_sei_av,
-            "X-averaged total "
-            + self.domain.lower()
-            + " sei thickness [m]": L_sei_av * L_scale,
-            "Inner "
-            + self.domain.lower()
-            + " sei concentration [mol.m-3]": n_inner * n_scale,
+            "Inner " + domain + " sei thickness": L_inner,
+            "Inner " + domain + " sei thickness [m]": L_inner * L_scale,
+            "X-averaged inner " + domain + " sei thickness": L_inner_av,
+            "X-averaged inner " + domain + " sei thickness [m]": L_inner_av * L_scale,
+            "Outer " + domain + " sei thickness": L_outer,
+            "Outer " + domain + " sei thickness [m]": L_outer * L_scale,
+            "X-averaged outer " + domain + " sei thickness": L_outer_av,
+            "X-averaged outer " + domain + " sei thickness [m]": L_outer_av * L_scale,
+            "Total " + domain + " sei thickness": L_sei,
+            "Total " + domain + " sei thickness [m]": L_sei * L_scale,
+            "X-averaged total " + domain + " sei thickness": L_sei_av,
+            "X-averaged total " + domain + " sei thickness [m]": L_sei_av * L_scale,
+            "Inner " + domain + " sei concentration [mol.m-3]": n_inner * n_scale,
             "X-averaged inner "
-            + self.domain.lower()
+            + domain
             + " sei concentration [mol.m-3]": n_inner_av * n_scale,
-            "Outer "
-            + self.domain.lower()
-            + " sei concentration [mol.m-3]": n_outer * n_outer_scale,
+            "Outer " + domain + " sei concentration [mol.m-3]": n_outer * n_outer_scale,
             "X-averaged outer "
-            + self.domain.lower()
+            + domain
             + " sei concentration [mol.m-3]": n_outer_av * n_outer_scale,
             self.domain + " sei concentration [mol.m-3]": n_SEI * n_scale,
-            "X-averaged "
-            + self.domain.lower()
-            + " sei concentration [mol.m-3]": n_SEI_av * n_scale,
-            "Loss of lithium to "
-            + self.domain.lower()
-            + " sei [mols]": Q_sei * n_scale,
+            "X-averaged " + domain + " sei concentration [mol.m-3]": n_SEI_av * n_scale,
+            "Loss of lithium to " + domain + " sei [mols]": Q_sei * n_scale,
         }
 
         return variables
@@ -145,41 +134,34 @@ def _get_standard_reaction_variables(self, j_inner, j_outer):
         j_sei = j_inner + j_outer
         j_sei_av = pybamm.x_average(j_sei)
 
+        domain = self.domain.lower() + " electrode"
+        Domain = domain.capitalize()
+
         variables = {
+            "Inner " + domain + " sei interfacial current density": j_inner,
             "Inner "
-            + self.domain.lower()
-            + " sei interfacial current density": j_inner,
-            "Inner "
-            + self.domain.lower()
+            + domain
             + " sei interfacial current density [A.m-2]": j_inner * j_scale,
+            "X-averaged inner " + domain + " sei interfacial current density": j_i_av,
             "X-averaged inner "
-            + self.domain.lower()
-            + " sei interfacial current density": j_i_av,
-            "X-averaged inner "
-            + self.domain.lower()
+            + domain
             + " sei interfacial current density [A.m-2]": j_i_av * j_scale,
+            "Outer " + domain + " sei interfacial current density": j_outer,
             "Outer "
-            + self.domain.lower()
-            + " sei interfacial current density": j_outer,
-            "Outer "
-            + self.domain.lower()
+            + domain
             + " sei interfacial current density [A.m-2]": j_outer * j_scale,
+            "X-averaged outer " + domain + " sei interfacial current density": j_o_av,
             "X-averaged outer "
-            + self.domain.lower()
-            + " sei interfacial current density": j_o_av,
-            "X-averaged outer "
-            + self.domain.lower()
+            + domain
             + " sei interfacial current density [A.m-2]": j_o_av * j_scale,
-            self.domain + " sei interfacial current density": j_sei,
-            self.domain + " sei interfacial current density [A.m-2]": j_sei * j_scale,
-            "X-averaged "
-            + self.domain.lower()
-            + " sei interfacial current density": j_sei_av,
+            Domain + " sei interfacial current density": j_sei,
+            Domain + " sei interfacial current density [A.m-2]": j_sei * j_scale,
+            "X-averaged " + domain + " sei interfacial current density": j_sei_av,
             "X-averaged "
-            + self.domain.lower()
+            + domain
             + " sei interfacial current density [A.m-2]": j_sei_av * j_scale,
             "Scaled "
-            + self.domain.lower()
+            + domain
             + " sei interfacial current density": j_sei * Gamma_SEI_n,
         }
 
diff --git a/pybamm/models/submodels/interface/sei/constant_sei.py b/pybamm/models/submodels/interface/sei/constant_sei.py
new file mode 100644
index 0000000000..61af717dd4
--- /dev/null
+++ b/pybamm/models/submodels/interface/sei/constant_sei.py
@@ -0,0 +1,52 @@
+#
+# Class for constant SEI thickness
+#
+import pybamm
+from .base_sei import BaseModel
+
+
+class ConstantSEI(BaseModel):
+    """Base class for SEI with constant thickness.
+
+    Note that there is no SEI current, so we don't need to update the "sum of
+    interfacial current densities" variables from
+    :class:`pybamm.interface.BaseInterface`
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    domain : str
+        The domain of the model either 'Negative' or 'Positive'
+
+    **Extends:** :class:`pybamm.sei.BaseModel`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def get_fundamental_variables(self):
+        # Constant thicknesses
+        L_inner = pybamm.sei_parameters.L_inner_0
+        L_outer = pybamm.sei_parameters.L_outer_0
+        variables = self._get_standard_thickness_variables(L_inner, L_outer)
+
+        # Reactions
+        zero = pybamm.FullBroadcast(
+            pybamm.Scalar(0), self.domain.lower() + " electrode", "current collector"
+        )
+        variables.update(self._get_standard_reaction_variables(zero, zero))
+
+        return variables
+
+    def get_coupled_variables(self, variables):
+        # Update whole cell variables, which also updates the "sum of" variables
+        if (
+            "Negative electrode sei interfacial current density" in variables
+            and "Positive electrode sei interfacial current density" in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+
+        return variables
diff --git a/pybamm/models/submodels/sei/electron_migration_limited.py b/pybamm/models/submodels/interface/sei/electron_migration_limited.py
similarity index 57%
rename from pybamm/models/submodels/sei/electron_migration_limited.py
rename to pybamm/models/submodels/interface/sei/electron_migration_limited.py
index 2ba791ec62..fbc6ba9467 100644
--- a/pybamm/models/submodels/sei/electron_migration_limited.py
+++ b/pybamm/models/submodels/interface/sei/electron_migration_limited.py
@@ -15,8 +15,7 @@ class ElectronMigrationLimited(BaseModel):
     domain : str
         The domain of the model either 'Negative' or 'Positive'
 
-
-    **Extends:** :class:`pybamm.particle.BaseParticle`
+    **Extends:** :class:`pybamm.sei.BaseModel`
     """
 
     def __init__(self, param, domain):
@@ -31,8 +30,10 @@ def get_fundamental_variables(self):
         return variables
 
     def get_coupled_variables(self, variables):
-        L_sei_inner = variables["Inner negative electrode sei thickness"]
-        phi_s_n = variables["Negative electrode potential"]
+        L_sei_inner = variables[
+            "Inner " + self.domain.lower() + " electrode sei thickness"
+        ]
+        phi_s_n = variables[self.domain + " electrode potential"]
 
         C_sei = pybamm.sei_parameters.C_sei_electron
         U_inner = pybamm.sei_parameters.U_inner_electron
@@ -45,25 +46,32 @@ def get_coupled_variables(self, variables):
 
         variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
 
+        # Update whole cell variables, which also updates the "sum of" variables
+        if (
+            "Negative electrode sei interfacial current density" in variables
+            and "Positive electrode sei interfacial current density" in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+
         return variables
 
     def set_rhs(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
-        j_inner = variables[
-            "Inner " + self.domain.lower() + " sei interfacial current density"
-        ]
-        j_outer = variables[
-            "Outer " + self.domain.lower() + " sei interfacial current density"
-        ]
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
+        j_inner = variables["Inner " + domain + " sei interfacial current density"]
+        j_outer = variables["Outer " + domain + " sei interfacial current density"]
 
         v_bar = pybamm.sei_parameters.v_bar
 
         self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
 
     def set_initial_conditions(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
 
         L_inner_0 = pybamm.sei_parameters.L_inner_0
         L_outer_0 = pybamm.sei_parameters.L_outer_0
diff --git a/pybamm/models/submodels/sei/interstitial_diffusion_limited.py b/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
similarity index 55%
rename from pybamm/models/submodels/sei/interstitial_diffusion_limited.py
rename to pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
index 4e78c13c91..0b2b1e68dc 100644
--- a/pybamm/models/submodels/sei/interstitial_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
@@ -15,8 +15,7 @@ class InterstitialDiffusionLimited(BaseModel):
     domain : str
         The domain of the model either 'Negative' or 'Positive'
 
-
-    **Extends:** :class:`pybamm.particle.BaseParticle`
+    **Extends:** :class:`pybamm.sei.BaseModel`
     """
 
     def __init__(self, param, domain):
@@ -31,9 +30,11 @@ def get_fundamental_variables(self):
         return variables
 
     def get_coupled_variables(self, variables):
-        L_sei_inner = variables["Inner negative electrode sei thickness"]
-        phi_s_n = variables["Negative electrode potential"]
-        phi_e_n = variables["Negative electrolyte potential"]
+        L_sei_inner = variables[
+            "Inner " + self.domain.lower() + " electrode sei thickness"
+        ]
+        phi_s_n = variables[self.domain + " electrode potential"]
+        phi_e_n = variables[self.domain + " electrolyte potential"]
 
         C_sei = pybamm.sei_parameters.C_sei_inter
 
@@ -45,25 +46,32 @@ def get_coupled_variables(self, variables):
 
         variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
 
+        # Update whole cell variables, which also updates the "sum of" variables
+        if (
+            "Negative electrode sei interfacial current density" in variables
+            and "Positive electrode sei interfacial current density" in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+
         return variables
 
     def set_rhs(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
-        j_inner = variables[
-            "Inner " + self.domain.lower() + " sei interfacial current density"
-        ]
-        j_outer = variables[
-            "Outer " + self.domain.lower() + " sei interfacial current density"
-        ]
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
+        j_inner = variables["Inner " + domain + " sei interfacial current density"]
+        j_outer = variables["Outer " + domain + " sei interfacial current density"]
 
         v_bar = pybamm.sei_parameters.v_bar
 
         self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
 
     def set_initial_conditions(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
 
         L_inner_0 = pybamm.sei_parameters.L_inner_0
         L_outer_0 = pybamm.sei_parameters.L_outer_0
diff --git a/pybamm/models/submodels/sei/no_sei.py b/pybamm/models/submodels/interface/sei/no_sei.py
similarity index 52%
rename from pybamm/models/submodels/sei/no_sei.py
rename to pybamm/models/submodels/interface/sei/no_sei.py
index 65098b8713..17d659dacc 100644
--- a/pybamm/models/submodels/sei/no_sei.py
+++ b/pybamm/models/submodels/interface/sei/no_sei.py
@@ -15,8 +15,7 @@ class NoSEI(BaseModel):
     domain : str
         The domain of the model either 'Negative' or 'Positive'
 
-
-    **Extends:** :class:`pybamm.particle.BaseParticle`
+    **Extends:** :class:`pybamm.sei.BaseModel`
     """
 
     def __init__(self, param, domain):
@@ -24,8 +23,20 @@ def __init__(self, param, domain):
 
     def get_fundamental_variables(self):
         zero = pybamm.FullBroadcast(
-            pybamm.Scalar(0), self.domain.lower(), "current collector"
+            pybamm.Scalar(0), self.domain.lower() + " electrode", "current collector"
         )
         variables = self._get_standard_thickness_variables(zero, zero)
         variables.update(self._get_standard_reaction_variables(zero, zero))
         return variables
+
+    def get_coupled_variables(self, variables):
+        # Update whole cell variables, which also updates the "sum of" variables
+        if (
+            "Negative electrode sei interfacial current density" in variables
+            and "Positive electrode sei interfacial current density" in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+
+        return variables
diff --git a/pybamm/models/submodels/interface/sei/reaction_limited.py b/pybamm/models/submodels/interface/sei/reaction_limited.py
new file mode 100644
index 0000000000..b8a43d7bad
--- /dev/null
+++ b/pybamm/models/submodels/interface/sei/reaction_limited.py
@@ -0,0 +1,94 @@
+#
+# Class for reaction limited SEI growth
+#
+import pybamm
+from .base_sei import BaseModel
+
+
+class ReactionLimited(BaseModel):
+    """Base class for reaction limited SEI growth.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    domain : str
+        The domain of the model either 'Negative' or 'Positive'
+
+    **Extends:** :class:`pybamm.sei.BaseModel`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def get_fundamental_variables(self):
+        L_inner = pybamm.standard_variables.L_inner
+        L_outer = pybamm.standard_variables.L_outer
+
+        variables = self._get_standard_thickness_variables(L_inner, L_outer)
+
+        return variables
+
+    def get_coupled_variables(self, variables):
+        phi_s_n = variables[self.domain + " electrode potential"]
+        phi_e_n = variables[self.domain + " electrolyte potential"]
+
+        # Look for current that contributes to the -IR drop
+        # If we can't find the interfacial current density from the main reaction, j,
+        # it's ok to fall back on the total interfacial current density, j_tot
+        # This should only happen when the interface submodel is "InverseButlerVolmer"
+        # in which case j = j_tot (uniform) anyway
+        try:
+            j = variables[self.domain + " electrode interfacial current density"]
+        except KeyError:
+            j = variables[
+                "X-averaged "
+                + self.domain.lower()
+                + " electrode total interfacial current density"
+            ]
+        L_sei = variables["Total " + self.domain.lower() + " electrode sei thickness"]
+
+        C_sei = pybamm.sei_parameters.C_sei_reaction
+        R_sei = pybamm.sei_parameters.R_sei
+        alpha = 0.5
+        # alpha = pybamm.sei_parameters.alpha
+
+        # need to revise for thermal case
+        j_sei = -(1 / C_sei) * pybamm.exp(-(phi_s_n - phi_e_n - j * L_sei * R_sei))
+
+        j_inner = alpha * j_sei
+        j_outer = (1 - alpha) * j_sei
+
+        variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
+
+        # Update whole cell variables, which also updates the "sum of" variables
+        if (
+            "Negative electrode sei interfacial current density" in variables
+            and "Positive electrode sei interfacial current density" in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+
+        return variables
+
+    def set_rhs(self, variables):
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
+        j_inner = variables["Inner " + domain + " sei interfacial current density"]
+        j_outer = variables["Outer " + domain + " sei interfacial current density"]
+
+        v_bar = pybamm.sei_parameters.v_bar
+
+        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
+
+    def set_initial_conditions(self, variables):
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
+
+        L_inner_0 = pybamm.sei_parameters.L_inner_0
+        L_outer_0 = pybamm.sei_parameters.L_outer_0
+
+        self.initial_conditions = {L_inner: L_inner_0, L_outer: L_outer_0}
diff --git a/pybamm/models/submodels/sei/solvent_diffusion_limited.py b/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
similarity index 57%
rename from pybamm/models/submodels/sei/solvent_diffusion_limited.py
rename to pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
index 9b1b022421..6d7b3bd8cf 100644
--- a/pybamm/models/submodels/sei/solvent_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
@@ -15,8 +15,7 @@ class SolventDiffusionLimited(BaseModel):
     domain : str
         The domain of the model either 'Negative' or 'Positive'
 
-
-    **Extends:** :class:`pybamm.particle.BaseParticle`
+    **Extends:** :class:`pybamm.sei.BaseModel`
     """
 
     def __init__(self, param, domain):
@@ -31,7 +30,9 @@ def get_fundamental_variables(self):
         return variables
 
     def get_coupled_variables(self, variables):
-        L_sei_outer = variables["Outer negative electrode sei thickness"]
+        L_sei_outer = variables[
+            "Outer " + self.domain.lower() + " electrode sei thickness"
+        ]
 
         C_sei = pybamm.sei_parameters.C_sei_solvent
 
@@ -43,25 +44,32 @@ def get_coupled_variables(self, variables):
 
         variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
 
+        # Update whole cell variables, which also updates the "sum of" variables
+        if (
+            "Negative electrode sei interfacial current density" in variables
+            and "Positive electrode sei interfacial current density" in variables
+        ):
+            variables.update(
+                self._get_standard_whole_cell_interfacial_current_variables(variables)
+            )
+
         return variables
 
     def set_rhs(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
-        j_inner = variables[
-            "Inner " + self.domain.lower() + " sei interfacial current density"
-        ]
-        j_outer = variables[
-            "Outer " + self.domain.lower() + " sei interfacial current density"
-        ]
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
+        j_inner = variables["Inner " + domain + " sei interfacial current density"]
+        j_outer = variables["Outer " + domain + " sei interfacial current density"]
 
         v_bar = pybamm.sei_parameters.v_bar
 
         self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
 
     def set_initial_conditions(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
+        domain = self.domain.lower() + " electrode"
+        L_inner = variables["Inner " + domain + " sei thickness"]
+        L_outer = variables["Outer " + domain + " sei thickness"]
 
         L_inner_0 = pybamm.sei_parameters.L_inner_0
         L_outer_0 = pybamm.sei_parameters.L_outer_0
diff --git a/pybamm/models/submodels/sei/constant_sei.py b/pybamm/models/submodels/sei/constant_sei.py
deleted file mode 100644
index 43b450bee0..0000000000
--- a/pybamm/models/submodels/sei/constant_sei.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#
-# Class for constant SEI thickness
-#
-import pybamm
-from .base_sei import BaseModel
-
-
-class ConstantSEI(BaseModel):
-    """Base class for no SEI.
-
-    Parameters
-    ----------
-    param : parameter class
-        The parameters to use for this submodel
-    domain : str
-        The domain of the model either 'Negative' or 'Positive'
-
-    **Extends:** :class:`pybamm.particle.BaseParticle`
-    """
-
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
-
-    def get_fundamental_variables(self):
-        param = self.param
-
-        # Constant thicknesses
-        L_inner = pybamm.sei_parameters.L_inner_0
-        L_outer = pybamm.sei_parameters.L_outer_0
-        variables = self._get_standard_thickness_variables(L_inner, L_outer)
-
-        # Reactions
-        zero = pybamm.FullBroadcast(
-            pybamm.Scalar(0), self.domain.lower(), "current collector"
-        )
-        variables.update(self._get_standard_reaction_variables(zero, zero))
-
-        return variables
diff --git a/pybamm/models/submodels/sei/reaction_limited.py b/pybamm/models/submodels/sei/reaction_limited.py
deleted file mode 100644
index 4fe2f0f6ad..0000000000
--- a/pybamm/models/submodels/sei/reaction_limited.py
+++ /dev/null
@@ -1,75 +0,0 @@
-#
-# Class for reaction limited SEI growth
-#
-import pybamm
-from .base_sei import BaseModel
-
-
-class ReactionLimited(BaseModel):
-    """Base class for reaction limited SEI growth.
-
-    Parameters
-    ----------
-    param : parameter class
-        The parameters to use for this submodel
-    domain : str
-        The domain of the model either 'Negative' or 'Positive'
-
-
-    **Extends:** :class:`pybamm.particle.BaseParticle`
-    """
-
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
-
-    def get_fundamental_variables(self):
-        L_inner = pybamm.standard_variables.L_inner
-        L_outer = pybamm.standard_variables.L_outer
-
-        variables = self._get_standard_thickness_variables(L_inner, L_outer)
-
-        return variables
-
-    def get_coupled_variables(self, variables):
-        phi_s_n = variables["Negative electrode potential"]
-        phi_e_n = variables["Negative electrolyte potential"]
-        j_n = variables["Negative electrode interfacial current density"]
-        L_sei = variables["Total negative electrode sei thickness"]
-
-        C_sei = pybamm.sei_parameters.C_sei_reaction
-        R_sei = pybamm.sei_parameters.R_sei
-        alpha = 0.5
-        # alpha = pybamm.sei_parameters.alpha
-
-        # need to revise for thermal case
-        j_sei = - (1 / C_sei) * pybamm.exp(-(phi_s_n - phi_e_n - j_n * L_sei * R_sei))
-
-        j_inner = alpha * j_sei
-        j_outer = (1 - alpha) * j_sei
-
-        variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
-
-        return variables
-
-    def set_rhs(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
-        j_inner = variables[
-            "Inner " + self.domain.lower() + " sei interfacial current density"
-        ]
-        j_outer = variables[
-            "Outer " + self.domain.lower() + " sei interfacial current density"
-        ]
-
-        v_bar = pybamm.sei_parameters.v_bar
-
-        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
-
-    def set_initial_conditions(self, variables):
-        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
-        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
-
-        L_inner_0 = pybamm.sei_parameters.L_inner_0
-        L_outer_0 = pybamm.sei_parameters.L_outer_0
-
-        self.initial_conditions = {L_inner: L_inner_0, L_outer: L_outer_0}
diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py
index c4f830c8bf..2f7608ea4e 100644
--- a/pybamm/parameters/standard_parameters_lithium_ion.py
+++ b/pybamm/parameters/standard_parameters_lithium_ion.py
@@ -97,8 +97,6 @@
 C_dl_p_dimensional = pybamm.Parameter(
     "Positive electrode double-layer capacity [F.m-2]"
 )
-# Oxygen parameters, for reusing same submodels as lead-acid
-s_plus_Ox = 0
 
 
 # Initial conditions
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index 8a04484ca6..f148e72aad 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -107,7 +107,7 @@ def test_default_spatial_methods(self):
             )
         )
 
-    def test_bad_options(self):
+    def test_options(self):
         with self.assertRaisesRegex(pybamm.OptionError, "Option"):
             pybamm.BaseBatteryModel({"bad option": "bad option"})
         with self.assertRaisesRegex(pybamm.OptionError, "current collector model"):
@@ -131,6 +131,17 @@ def test_bad_options(self):
         with self.assertRaisesRegex(pybamm.OptionError, "operating mode"):
             pybamm.BaseBatteryModel({"operating mode": "bad operating mode"})
 
+        # SEI options
+        with self.assertRaisesRegex(pybamm.OptionError, "sei"):
+            pybamm.BaseBatteryModel({"sei": "bad sei"})
+        with self.assertRaisesRegex(pybamm.OptionError, "sei film resistance"):
+            pybamm.BaseBatteryModel({"sei film resistance": "bad sei film resistance"})
+        # variable defaults
+        model = pybamm.BaseBatteryModel()
+        self.assertEqual(model.options["sei film resistance"], None)
+        model = pybamm.BaseBatteryModel({"sei": "constant"})
+        self.assertEqual(model.options["sei film resistance"], "distributed")
+
     def test_build_twice(self):
         model = pybamm.lithium_ion.SPM()  # need to pick a model to set vars and build
         with self.assertRaisesRegex(pybamm.ModelError, "Model already built"):
@@ -141,7 +152,7 @@ def test_get_coupled_variables(self):
         model.submodels["current collector"] = pybamm.current_collector.Uniform(
             model.param
         )
-        with self.assertRaisesRegex(pybamm.ModelError, "Submodel"):
+        with self.assertRaisesRegex(pybamm.ModelError, "Missing variable"):
             model.build_model()
 
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_base_lead_acid_model.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_base_lead_acid_model.py
index dc0d7b7596..fdd0e859d2 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_base_lead_acid_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_base_lead_acid_model.py
@@ -30,9 +30,11 @@ def test_default_geometry(self):
     def test_incompatible_options(self):
         with self.assertRaisesRegex(
             pybamm.OptionError,
-            "Lead-acid models can only have thermal " "effects if dimensionality is 0.",
+            "Lead-acid models can only have thermal effects if dimensionality is 0.",
         ):
             pybamm.lead_acid.BaseModel({"dimensionality": 1, "thermal": "x-full"})
+        with self.assertRaisesRegex(pybamm.OptionError, "SEI"):
+            pybamm.lead_acid.BaseModel({"sei": "constant"})
 
 
 if __name__ == "__main__":

From 1f74cddbdcd0e59af676daccb4f5036b703aa330 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 13 May 2020 16:35:30 -0400
Subject: [PATCH 03/15] #984 add comments

---
 examples/scripts/sei_growth.py                | 21 ++++++++++++-------
 .../inverse_kinetics/inverse_butler_volmer.py | 19 ++++++++++++++++-
 2 files changed, 31 insertions(+), 9 deletions(-)

diff --git a/examples/scripts/sei_growth.py b/examples/scripts/sei_growth.py
index c7c9c4f088..7d622b6e4a 100644
--- a/examples/scripts/sei_growth.py
+++ b/examples/scripts/sei_growth.py
@@ -1,14 +1,19 @@
 import pybamm as pb
 import numpy as np
 
-pb.set_logging_level("DEBUG")
-
-options = {
-    "sei": "reaction limited",
-    "sei film resistance": None,
-    # "surface form": "algebraic",
-}
-models = [pb.lithium_ion.SPM(options), pb.lithium_ion.DFN(options)]
+pb.set_logging_level("INFO")
+
+models = [
+    pb.lithium_ion.SPM({"sei": "reaction limited", "sei film resistance": None}),
+    pb.lithium_ion.SPM(
+        {
+            "sei": "reaction limited",
+            "sei film resistance": None,
+            "surface form": "algebraic",
+        }
+    ),
+    pb.lithium_ion.DFN({"sei": "reaction limited", "sei film resistance": None}),
+]
 
 sims = []
 for model in models:
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 da59aa9740..59698ff531 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -59,6 +59,9 @@ def get_coupled_variables(self, variables):
         # eta_r is the overpotential from inverting Butler-Volmer, regardless of any
         # additional SEI resistance. What changes is how delta_phi is defined in terms
         # of eta_r
+        # We use the total resistance to calculate eta_r, but this only introduces
+        # negligible errors. For the exact answer, the surface form submodels should
+        # be used instead
         eta_r = self._get_overpotential(j_tot, j0, ne, T)
 
         # With SEI resistance (distributed and averaged have the same effect here)
@@ -89,7 +92,21 @@ def _get_overpotential(self, j, j0, ne, T):
 
 class CurrentForInverseButlerVolmer(BaseInterface):
     """
-    Submodel for the current associated with the inverse Butler-Volmer formulation. 
+    Submodel for the current associated with the inverse Butler-Volmer formulation. This
+    has to be created as a separate submodel because of how the interfacial currents
+    are calculated:
+
+    1. Calculate eta_r from the total average current j_tot_av = I_app / L
+    2. Calculate j_sei from eta_r
+    3. Calculate j = j_tot_av - j_sei
+
+    To be able to do step 1, then step 2, then step 3 requires two different submodels
+    for step 1 and step 2
+
+    This introduces a little bit of error because eta_r is calculated using j_tot_av
+    instead of j. But since j_sei is very small, this error is very small. The "surface
+    form" model solves a differential or algebraic equation for delta_phi, which gives
+    the exact right answer. Comparing the two approaches shows almost no difference.
 
     Parameters
     ----------

From 28fc8f1eacc69010fdeb473788d231d90b38b0ba Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 13 May 2020 18:11:23 -0400
Subject: [PATCH 04/15] #984 passing unit tests

---
 examples/scripts/compare_lithium_ion.py       | 12 ++----
 examples/scripts/sei_growth.py                |  9 ++--
 .../lead_acid/base_lead_acid_model.py         |  5 +++
 .../full_battery_models/lead_acid/full.py     |  1 +
 .../lead_acid/higher_order.py                 |  1 +
 .../full_battery_models/lead_acid/loqs.py     | 11 +++++
 .../full_battery_models/lithium_ion/spme.py   | 10 +++++
 .../base_electrolyte_diffusion.py             |  3 +-
 .../submodels/interface/base_interface.py     | 34 ++++++---------
 .../inverse_kinetics/inverse_butler_volmer.py |  1 +
 .../interface/kinetics/base_kinetics.py       | 41 ++++++++++---------
 .../submodels/interface/sei/constant_sei.py   |  1 +
 .../sei/electron_migration_limited.py         |  1 +
 .../sei/interstitial_diffusion_limited.py     |  1 +
 .../models/submodels/interface/sei/no_sei.py  |  1 +
 .../interface/sei/reaction_limited.py         |  1 +
 .../sei/solvent_diffusion_limited.py          |  1 +
 .../standard_parameters_lithium_ion.py        |  7 ++--
 .../test_lithium_ion/test_dfn.py              |  6 +++
 .../test_lithium_ion/test_dfn.py              |  5 +++
 20 files changed, 94 insertions(+), 58 deletions(-)

diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py
index c92dfe2aa5..42cd480390 100644
--- a/examples/scripts/compare_lithium_ion.py
+++ b/examples/scripts/compare_lithium_ion.py
@@ -16,21 +16,17 @@
     pybamm.set_logging_level("INFO")
 
 # load models
-# options = {"thermal": "lumped"}
+options = {"thermal": "lumped"}
 models = [
-    # pybamm.lithium_ion.SPMe(),
-    # pybamm.lithium_ion.SPMe({"sei": "constant"}),
-    # pybamm.lithium_ion.SPMe(options),
-    pybamm.lithium_ion.DFN(),
-    pybamm.lithium_ion.DFN({"sei": "constant", "sei film resistance": "average"}),
-    pybamm.lithium_ion.DFN({"sei": "constant", "sei film resistance": "distributed"}),
+    pybamm.lithium_ion.SPM(options),
+    pybamm.lithium_ion.SPMe(options),
+    pybamm.lithium_ion.DFN(options),
 ]
 
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
 param["Current function [A]"] = 1
-param["SEI resistivity [Ohm.m]"] = 1
 
 for model in models:
     param.process_model(model)
diff --git a/examples/scripts/sei_growth.py b/examples/scripts/sei_growth.py
index 7d622b6e4a..0391bef490 100644
--- a/examples/scripts/sei_growth.py
+++ b/examples/scripts/sei_growth.py
@@ -4,15 +4,16 @@
 pb.set_logging_level("INFO")
 
 models = [
-    pb.lithium_ion.SPM({"sei": "reaction limited", "sei film resistance": None}),
+    pb.lithium_ion.SPM({"sei": "reaction limited"}),
     pb.lithium_ion.SPM(
         {
             "sei": "reaction limited",
-            "sei film resistance": None,
+            # "sei film resistance": "average",
             "surface form": "algebraic",
-        }
+        },
+        name="Algebraic SPM",
     ),
-    pb.lithium_ion.DFN({"sei": "reaction limited", "sei film resistance": None}),
+    pb.lithium_ion.DFN({"sei": "reaction limited"}),
 ]
 
 sims = []
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 990434414c..826374bd99 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
@@ -58,3 +58,8 @@ def set_soc_variables(self):
             self.variables["Fractional Charge Input"] = fci
             self.rhs[fci] = -self.variables["Total current density"] * 100
             self.initial_conditions[fci] = self.param.q_init * 100
+
+    def set_sei_submodel(self):
+
+        self.submodels["negative sei"] = pybamm.sei.NoSEI(self.param, "Negative")
+        self.submodels["positive sei"] = pybamm.sei.NoSEI(self.param, "Positive")
diff --git a/pybamm/models/full_battery_models/lead_acid/full.py b/pybamm/models/full_battery_models/lead_acid/full.py
index 1fa1e70b15..3e4716e9bc 100644
--- a/pybamm/models/full_battery_models/lead_acid/full.py
+++ b/pybamm/models/full_battery_models/lead_acid/full.py
@@ -44,6 +44,7 @@ def __init__(self, options=None, name="Full model", build=True):
         self.set_thermal_submodel()
         self.set_side_reaction_submodels()
         self.set_current_collector_submodel()
+        self.set_sei_submodel()
 
         if build:
             self.build_model()
diff --git a/pybamm/models/full_battery_models/lead_acid/higher_order.py b/pybamm/models/full_battery_models/lead_acid/higher_order.py
index 672250a4aa..24a0b02eb7 100644
--- a/pybamm/models/full_battery_models/lead_acid/higher_order.py
+++ b/pybamm/models/full_battery_models/lead_acid/higher_order.py
@@ -52,6 +52,7 @@ def __init__(self, options=None, name="Composite model", build=True):
         self.set_tortuosity_submodels()
         self.set_thermal_submodel()
         self.set_current_collector_submodel()
+        self.set_sei_submodel()
 
         if build:
             self.build_model()
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index e40b860754..97cafcb6e3 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -44,6 +44,7 @@ def __init__(self, options=None, name="LOQS model", build=True):
         self.set_thermal_submodel()
         self.set_side_reaction_submodels()
         self.set_current_collector_submodel()
+        self.set_sei_submodel()
 
         if build:
             self.build_model()
@@ -140,6 +141,16 @@ def set_interfacial_submodel(self):
             ] = pybamm.interface.InverseButlerVolmer(
                 self.param, "Positive", "lead-acid main"
             )
+            self.submodels[
+                "negative interface current"
+            ] = pybamm.interface.CurrentForInverseButlerVolmer(
+                self.param, "Negative", "lithium-ion main"
+            )
+            self.submodels[
+                "positive interface current"
+            ] = pybamm.interface.CurrentForInverseButlerVolmer(
+                self.param, "Positive", "lithium-ion main"
+            )
         else:
             self.submodels[
                 "leading-order negative interface"
diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py
index e033499988..93cd6deb44 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spme.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spme.py
@@ -83,6 +83,16 @@ def set_interfacial_submodel(self):
         self.submodels["positive interface"] = pybamm.interface.InverseButlerVolmer(
             self.param, "Positive", "lithium-ion main", self.options
         )
+        self.submodels[
+            "negative interface current"
+        ] = pybamm.interface.CurrentForInverseButlerVolmer(
+            self.param, "Negative", "lithium-ion main"
+        )
+        self.submodels[
+            "positive interface current"
+        ] = pybamm.interface.CurrentForInverseButlerVolmer(
+            self.param, "Positive", "lithium-ion main"
+        )
 
     def set_particle_submodel(self):
 
diff --git a/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py
index 42f016a281..7748cde679 100644
--- a/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py
+++ b/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py
@@ -96,8 +96,7 @@ def _get_standard_flux_variables(self, N_e):
         """
 
         param = self.param
-        D_e_typ = param.D_e(param.c_e_typ, param.T_init)
-        flux_scale = D_e_typ * param.c_e_typ / param.L_x
+        flux_scale = param.D_e_typ * param.c_e_typ / param.L_x
 
         variables = {
             "Electrolyte flux": N_e,
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index a6ed54dbd5..75f9346cdf 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -25,16 +25,19 @@ def __init__(self, param, domain, reaction):
         super().__init__(param, domain)
         if reaction == "lithium-ion main":
             self.reaction_name = ""  # empty reaction name for the main reaction
+            self.Reaction_icd = "Interfacial current density"
         elif reaction == "lead-acid main":
             self.reaction_name = ""  # empty reaction name for the main reaction
+            self.Reaction_icd = "Interfacial current density"
         elif reaction == "lead-acid oxygen":
             self.reaction_name = " oxygen"
+            self.Reaction_icd = "Oxygen interfacial current density"
         elif reaction == "lithium-ion oxygen":
             self.reaction_name = " oxygen"
+            self.Reaction_icd = "Oxygen interfacial current density"
         elif reaction == "sei":
             self.reaction_name = " sei"
-        else:
-            raise ValueError("Reaction name '{}' not recognized".format(reaction))
+            self.Reaction_icd = "Sei interfacial current density"
         self.reaction = reaction
 
     def _get_exchange_current_density(self, variables):
@@ -316,26 +319,13 @@ def _get_standard_whole_cell_interfacial_current_variables(self, variables):
         j = pybamm.Concatenation(j_n, j_s, j_p)
         j_dim = pybamm.Concatenation(j_n_scale * j_n, j_s, j_p_scale * j_p)
 
-        if self.reaction_name == "":
-            variables.update(
-                {
-                    "Interfacial current density": j,
-                    "Interfacial current density [A.m-2]": j_dim,
-                    "Interfacial current density per volume [A.m-3]": i_typ / L_x * j,
-                }
-            )
-        else:
-            reaction_name = self.reaction_name[1:].capitalize()
-            variables.update(
-                {
-                    reaction_name + " interfacial current density": j,
-                    reaction_name + " interfacial current density [A.m-2]": j_dim,
-                    reaction_name
-                    + " interfacial current density per volume [A.m-3]": i_typ
-                    / L_x
-                    * j,
-                }
-            )
+        variables.update(
+            {
+                self.Reaction_icd: j,
+                self.Reaction_icd + " [A.m-2]": j_dim,
+                self.Reaction_icd + " per volume [A.m-3]": i_typ / L_x * j,
+            }
+        )
 
         s_n, s_p = self._get_electrolyte_reaction_signed_stoichiometry()
         s = pybamm.Concatenation(
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 59698ff531..fe084c068d 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -143,6 +143,7 @@ def get_coupled_variables(self, variables):
             + self.reaction_name
             + " interfacial current density"
             in variables
+            and self.Reaction_icd not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index 1f67c28629..f6c79ee968 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -42,7 +42,12 @@ def get_fundamental_variables(self):
                 auxiliary_domains={"secondary": "current collector"},
             )
 
-            variables = self._get_standard_interfacial_current_variables(j)
+            variables = {
+                self.domain
+                + " electrode"
+                + self.reaction_name
+                + " interfacial current density variable": j
+            }
             return variables
         else:
             return {}
@@ -71,7 +76,9 @@ def get_coupled_variables(self, variables):
             L_sei = variables[
                 "Total " + self.domain.lower() + " electrode sei thickness"
             ]
-            j = variables[self.domain + " electrode interfacial current density"]
+            j = variables[
+                self.domain + " electrode interfacial current density variable"
+            ]
             eta_r -= j * L_sei * pybamm.sei_parameters.R_sei
         elif self.options["sei film resistance"] == "average":
             L_sei = variables[
@@ -91,15 +98,7 @@ def get_coupled_variables(self, variables):
         # found by solving an algebraic equation
         # (In the "distributed SEI resistance" model, we have already defined j)
         j = self._get_kinetics(j0, ne, eta_r, T)
-        if (
-            self.options["sei film resistance"] == "distributed"
-            and "main" in self.reaction
-        ):
-            variables.update(
-                {self.domain + " electrode" + self.reaction_name + " reaction term": j}
-            )
-        else:
-            variables.update(self._get_standard_interfacial_current_variables(j))
+        variables.update(self._get_standard_interfacial_current_variables(j))
 
         variables.update(
             self._get_standard_total_interfacial_current_variables(j_tot_av)
@@ -115,6 +114,7 @@ def get_coupled_variables(self, variables):
             + self.reaction_name
             + " interfacial current density"
             in variables
+            and self.Reaction_icd not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
@@ -130,17 +130,20 @@ def set_algebraic(self, variables):
             self.options["sei film resistance"] == "distributed"
             and "main" in self.reaction
         ):
+            j_var = variables[
+                self.domain
+                + " electrode"
+                + self.reaction_name
+                + " interfacial current density variable"
+            ]
             j = variables[
                 self.domain
                 + " electrode"
                 + self.reaction_name
                 + " interfacial current density"
             ]
-            reaction = variables[
-                self.domain + " electrode" + self.reaction_name + " reaction term"
-            ]
-            # Algebraic equation to set j equal to the reaction term
-            self.algebraic[j] = j - reaction
+            # Algebraic equation to set the variable j_var equal to the reaction term j
+            self.algebraic[j_var] = j_var - j
 
     def set_initial_conditions(self, variables):
         if (
@@ -148,18 +151,18 @@ def set_initial_conditions(self, variables):
             and "main" in self.reaction
         ):
             param = self.param
-            j = variables[
+            j_var = variables[
                 self.domain
                 + " electrode"
                 + self.reaction_name
-                + " interfacial current density"
+                + " interfacial current density variable"
             ]
             if self.domain == "Negative":
                 j_av_init = param.current_with_time / param.l_n
             elif self.domain == "Positive":
                 j_av_init = -param.current_with_time / param.l_p
 
-            self.initial_conditions[j] = j_av_init
+            self.initial_conditions[j_var] = j_av_init
 
     def _get_dj_dc(self, variables):
         """
diff --git a/pybamm/models/submodels/interface/sei/constant_sei.py b/pybamm/models/submodels/interface/sei/constant_sei.py
index 61af717dd4..e12a9ce43e 100644
--- a/pybamm/models/submodels/interface/sei/constant_sei.py
+++ b/pybamm/models/submodels/interface/sei/constant_sei.py
@@ -44,6 +44,7 @@ def get_coupled_variables(self, variables):
         if (
             "Negative electrode sei interfacial current density" in variables
             and "Positive electrode sei interfacial current density" in variables
+            and "Sei interfacial current density" not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/sei/electron_migration_limited.py b/pybamm/models/submodels/interface/sei/electron_migration_limited.py
index fbc6ba9467..9d76dd80e3 100644
--- a/pybamm/models/submodels/interface/sei/electron_migration_limited.py
+++ b/pybamm/models/submodels/interface/sei/electron_migration_limited.py
@@ -50,6 +50,7 @@ def get_coupled_variables(self, variables):
         if (
             "Negative electrode sei interfacial current density" in variables
             and "Positive electrode sei interfacial current density" in variables
+            and "Sei interfacial current density" not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py b/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
index 0b2b1e68dc..9e15caa522 100644
--- a/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
@@ -50,6 +50,7 @@ def get_coupled_variables(self, variables):
         if (
             "Negative electrode sei interfacial current density" in variables
             and "Positive electrode sei interfacial current density" in variables
+            and "Sei interfacial current density" not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/sei/no_sei.py b/pybamm/models/submodels/interface/sei/no_sei.py
index 17d659dacc..b23b951210 100644
--- a/pybamm/models/submodels/interface/sei/no_sei.py
+++ b/pybamm/models/submodels/interface/sei/no_sei.py
@@ -34,6 +34,7 @@ def get_coupled_variables(self, variables):
         if (
             "Negative electrode sei interfacial current density" in variables
             and "Positive electrode sei interfacial current density" in variables
+            and "Sei interfacial current density" not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/sei/reaction_limited.py b/pybamm/models/submodels/interface/sei/reaction_limited.py
index b8a43d7bad..89a250949b 100644
--- a/pybamm/models/submodels/interface/sei/reaction_limited.py
+++ b/pybamm/models/submodels/interface/sei/reaction_limited.py
@@ -65,6 +65,7 @@ def get_coupled_variables(self, variables):
         if (
             "Negative electrode sei interfacial current density" in variables
             and "Positive electrode sei interfacial current density" in variables
+            and "Sei interfacial current density" not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py b/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
index 6d7b3bd8cf..a1e638ba34 100644
--- a/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
@@ -48,6 +48,7 @@ def get_coupled_variables(self, variables):
         if (
             "Negative electrode sei interfacial current density" in variables
             and "Positive electrode sei interfacial current density" in variables
+            and "Sei interfacial current density" not in variables
         ):
             variables.update(
                 self._get_standard_whole_cell_interfacial_current_variables(variables)
diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py
index 2f7608ea4e..12dcd308e3 100644
--- a/pybamm/parameters/standard_parameters_lithium_ion.py
+++ b/pybamm/parameters/standard_parameters_lithium_ion.py
@@ -253,7 +253,8 @@ def U_p_dimensional(sto, T):
 tau_r_p = F * c_p_max / (j0_p_ref_dimensional * a_p_dim)
 
 # Electrolyte diffusion timescale
-tau_diffusion_e = L_x ** 2 / D_e_dimensional(c_e_typ, T_ref)
+D_e_typ = D_e_dimensional(c_e_typ, T_ref)
+tau_diffusion_e = L_x ** 2 / D_e_typ
 
 # Particle diffusion timescales
 tau_diffusion_n = R_n ** 2 / D_n_dimensional(pybamm.Scalar(1), T_ref)
@@ -434,13 +435,13 @@ def D_e(c_e, T):
     "Dimensionless electrolyte diffusivity"
     c_e_dimensional = c_e * c_e_typ
     T_dim = Delta_T * T + T_ref
-    return D_e_dimensional(c_e_dimensional, T_dim) / D_e_dimensional(c_e_typ, T_ref)
+    return D_e_dimensional(c_e_dimensional, T_dim) / D_e_typ
 
 
 def kappa_e(c_e, T):
     "Dimensionless electrolyte conductivity"
     c_e_dimensional = c_e * c_e_typ
-    kappa_scale = F ** 2 * D_e_dimensional(c_e_typ, T_ref) * c_e_typ / (R * T_ref)
+    kappa_scale = F ** 2 * D_e_typ * c_e_typ / (R * T_ref)
     T_dim = Delta_T * T + T_ref
     return kappa_e_dimensional(c_e_dimensional, T_dim) / kappa_scale
 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 2053f56798..8bfe950786 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -130,6 +130,12 @@ def test_well_posed_reaction_limited(self):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_well_posed_reaction_limited_average_film_resistance(self):
+        options = {"sei": "reaction limited", "sei film resistance": "average"}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
     def test_well_posed_solvent_diffusion_limited(self):
         options = {"sei": "solvent-diffusion limited"}
         model = pybamm.lithium_ion.DFN(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 869e489254..a4d201b8a8 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -124,6 +124,11 @@ def test_well_posed_reaction_limited(self):
         model = pybamm.lithium_ion.DFN(options)
         model.check_well_posedness()
 
+    def test_well_posed_reaction_limited_average_film_resistance(self):
+        options = {"sei": "reaction limited", "sei film resistance": "average"}
+        model = pybamm.lithium_ion.DFN(options)
+        model.check_well_posedness()
+
     def test_well_posed_solvent_diffusion_limited(self):
         options = {"sei": "solvent-diffusion limited"}
         model = pybamm.lithium_ion.DFN(options)

From 09c9f95f2b10b8ea45727aeed6bfbdb3c4387f1b Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 13 May 2020 21:50:52 -0400
Subject: [PATCH 05/15] #984 fix tests

---
 examples/notebooks/using-submodels.ipynb      | 173 +++++++++++-------
 .../{sei_growth.py => calendar_ageing.py}     |   7 +-
 examples/scripts/compare_lead_acid.py         |   4 +
 examples/scripts/custom_model.py              |  12 ++
 .../full_battery_models/lead_acid/loqs.py     |   4 +-
 .../inverse_kinetics/inverse_butler_volmer.py |   4 +-
 .../test_models/standard_output_tests.py      |  16 +-
 7 files changed, 138 insertions(+), 82 deletions(-)
 rename examples/scripts/{sei_growth.py => calendar_ageing.py} (90%)

diff --git a/examples/notebooks/using-submodels.ipynb b/examples/notebooks/using-submodels.ipynb
index ac6fc3d902..8bf2043c93 100644
--- a/examples/notebooks/using-submodels.ipynb
+++ b/examples/notebooks/using-submodels.ipynb
@@ -40,26 +40,7 @@
    "cell_type": "code",
    "execution_count": 2,
    "metadata": {},
-   "outputs": [
-    {
-     "ename": "DomainError",
-     "evalue": "Primary broadcast from current collector domain must be to electrode\n                or separator",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mDomainError\u001b[0m                               Traceback (most recent call last)",
-      "\u001b[0;32m<ipython-input-2-64212570b990>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpybamm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlithium_ion\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSPM\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/models/full_battery_models/lithium_ion/spm.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, options, name, build)\u001b[0m\n\u001b[1;32m     46\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mbuild\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 48\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     49\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mset_porosity_submodel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/models/full_battery_models/base_battery_model.py\u001b[0m in \u001b[0;36mbuild_model\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    493\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_fundamental_and_external\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    494\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 495\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_coupled_variables\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    496\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    497\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_model_equations\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/models/full_battery_models/base_battery_model.py\u001b[0m in \u001b[0;36mbuild_coupled_variables\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    425\u001b[0m                     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    426\u001b[0m                         self.variables.update(\n\u001b[0;32m--> 427\u001b[0;31m                             \u001b[0msubmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_coupled_variables\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    428\u001b[0m                         )\n\u001b[1;32m    429\u001b[0m                         \u001b[0msubmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mremove\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msubmodel_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/models/submodels/particle/fickian/fickian_single_particle.py\u001b[0m in \u001b[0;36mget_coupled_variables\u001b[0;34m(self, variables)\u001b[0m\n\u001b[1;32m     44\u001b[0m         T_k_xav = pybamm.PrimaryBroadcast(\n\u001b[1;32m     45\u001b[0m             \u001b[0mvariables\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"X-averaged \"\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdomain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlower\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m\" electrode temperature\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 46\u001b[0;31m             \u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdomain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlower\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m\" particle\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     47\u001b[0m         )\n\u001b[1;32m     48\u001b[0m         \u001b[0mN_s_xav\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flux_law\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc_s_xav\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mT_k_xav\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/expression_tree/broadcasts.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, child, broadcast_domain, name)\u001b[0m\n\u001b[1;32m     85\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     86\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchild\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbroadcast_domain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 87\u001b[0;31m         \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mchild\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbroadcast_domain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbroadcast_type\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"primary\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     88\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m     def check_and_set_domains(\n",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/expression_tree/broadcasts.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, child, broadcast_domain, broadcast_auxiliary_domains, broadcast_type, name)\u001b[0m\n\u001b[1;32m     53\u001b[0m         \u001b[0;31m# perform some basic checks and set attributes\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     54\u001b[0m         domain, auxiliary_domains = self.check_and_set_domains(\n\u001b[0;32m---> 55\u001b[0;31m             \u001b[0mchild\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbroadcast_type\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbroadcast_domain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbroadcast_auxiliary_domains\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     56\u001b[0m         )\n\u001b[1;32m     57\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbroadcast_type\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbroadcast_type\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
-      "\u001b[0;32m~/Documents/Energy_storage/PyBaMM/pybamm/expression_tree/broadcasts.py\u001b[0m in \u001b[0;36mcheck_and_set_domains\u001b[0;34m(self, child, broadcast_type, broadcast_domain, broadcast_auxiliary_domains)\u001b[0m\n\u001b[1;32m    102\u001b[0m             raise pybamm.DomainError(\n\u001b[1;32m    103\u001b[0m                 \"\"\"Primary broadcast from current collector domain must be to electrode\n\u001b[0;32m--> 104\u001b[0;31m                 or separator\"\"\"\n\u001b[0m\u001b[1;32m    105\u001b[0m             )\n\u001b[1;32m    106\u001b[0m         elif child.domain[0] in [\n",
-      "\u001b[0;31mDomainError\u001b[0m: Primary broadcast from current collector domain must be to electrode\n                or separator"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "model = pybamm.lithium_ion.SPM()"
    ]
@@ -73,28 +54,35 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "external circuit <pybamm.models.submodels.external_circuit.current_control_external_circuit.CurrentControl object at 0x13674b310>\n",
-      "porosity <pybamm.models.submodels.porosity.constant_porosity.Constant object at 0x13615a850>\n",
-      "electrolyte tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x13674b610>\n",
-      "electrode tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x13615a250>\n",
-      "convection <pybamm.models.submodels.convection.no_convection.NoConvection object at 0x1367f2e50>\n",
-      "negative interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x1367f2c90>\n",
-      "positive interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x1355c46d0>\n",
-      "negative particle <pybamm.models.submodels.particle.fickian_single_particle.FickianSingleParticle object at 0x1367e71d0>\n",
-      "positive particle <pybamm.models.submodels.particle.fickian_single_particle.FickianSingleParticle object at 0x1367f2ed0>\n",
-      "negative electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x1367f2e90>\n",
-      "leading-order electrolyte conductivity <pybamm.models.submodels.electrolyte_conductivity.leading_order_conductivity.LeadingOrder object at 0x1363dd590>\n",
-      "electrolyte diffusion <pybamm.models.submodels.electrolyte_diffusion.constant_concentration.ConstantConcentration object at 0x1367f2e10>\n",
-      "positive electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x1360b9890>\n",
-      "thermal <pybamm.models.submodels.thermal.isothermal.isothermal.Isothermal object at 0x1367f2f10>\n",
-      "current collector <pybamm.models.submodels.current_collector.homogeneous_current_collector.Uniform object at 0x136336e10>\n"
+      "external circuit <pybamm.models.submodels.external_circuit.current_control_external_circuit.CurrentControl object at 0x13e3106d0>\n",
+      "porosity <pybamm.models.submodels.porosity.constant_porosity.Constant object at 0x13e327110>\n",
+      "electrolyte tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x13e310510>\n",
+      "electrode tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x13e327390>\n",
+      "through-cell convection <pybamm.models.submodels.convection.through_cell.no_convection.NoConvection object at 0x13e327490>\n",
+      "transverse convection <pybamm.models.submodels.convection.transverse.no_convection.NoConvection object at 0x13e27c090>\n",
+      "negative interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x13e327650>\n",
+      "positive interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x13e327710>\n",
+      "negative interface current <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x13e327810>\n",
+      "positive interface current <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x13dc80950>\n",
+      "negative oxygen interface <pybamm.models.submodels.interface.kinetics.no_reaction.NoReaction object at 0x13e327950>\n",
+      "positive oxygen interface <pybamm.models.submodels.interface.kinetics.no_reaction.NoReaction object at 0x111501390>\n",
+      "negative particle <pybamm.models.submodels.particle.fickian_single_particle.FickianSingleParticle object at 0x111504050>\n",
+      "positive particle <pybamm.models.submodels.particle.fickian_single_particle.FickianSingleParticle object at 0x1114f5410>\n",
+      "negative electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x13df0ac10>\n",
+      "leading-order electrolyte conductivity <pybamm.models.submodels.electrolyte_conductivity.leading_order_conductivity.LeadingOrder object at 0x13dc80190>\n",
+      "electrolyte diffusion <pybamm.models.submodels.electrolyte_diffusion.constant_concentration.ConstantConcentration object at 0x13e327a10>\n",
+      "positive electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x13e327b10>\n",
+      "thermal <pybamm.models.submodels.thermal.isothermal.Isothermal object at 0x13e327c10>\n",
+      "current collector <pybamm.models.submodels.current_collector.homogeneous_current_collector.Uniform object at 0x13e327d10>\n",
+      "negative sei <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x13e327ed0>\n",
+      "positive sei <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x13e09d9d0>\n"
      ]
     }
    ],
@@ -112,7 +100,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -128,7 +116,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -144,28 +132,35 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "external circuit <pybamm.models.submodels.external_circuit.current_control_external_circuit.CurrentControl object at 0x136d34190>\n",
-      "porosity <pybamm.models.submodels.porosity.constant_porosity.Constant object at 0x136d34250>\n",
-      "electrolyte tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x136d34310>\n",
-      "electrode tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x136d343d0>\n",
-      "convection <pybamm.models.submodels.convection.no_convection.NoConvection object at 0x136d34410>\n",
-      "negative interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x136d34450>\n",
-      "positive interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x136d34490>\n",
-      "negative particle <pybamm.models.submodels.particle.fast_single_particle.FastSingleParticle object at 0x1367f2e10>\n",
-      "positive particle <pybamm.models.submodels.particle.fickian_single_particle.FickianSingleParticle object at 0x136d34510>\n",
-      "negative electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x136d34550>\n",
-      "leading-order electrolyte conductivity <pybamm.models.submodels.electrolyte_conductivity.leading_order_conductivity.LeadingOrder object at 0x136d34590>\n",
-      "electrolyte diffusion <pybamm.models.submodels.electrolyte_diffusion.constant_concentration.ConstantConcentration object at 0x136d345d0>\n",
-      "positive electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x136d34610>\n",
-      "thermal <pybamm.models.submodels.thermal.isothermal.isothermal.Isothermal object at 0x136d34650>\n",
-      "current collector <pybamm.models.submodels.current_collector.homogeneous_current_collector.Uniform object at 0x136d34690>\n"
+      "external circuit <pybamm.models.submodels.external_circuit.current_control_external_circuit.CurrentControl object at 0x13e80c550>\n",
+      "porosity <pybamm.models.submodels.porosity.constant_porosity.Constant object at 0x13e80c710>\n",
+      "electrolyte tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x13e80c990>\n",
+      "electrode tortuosity <pybamm.models.submodels.tortuosity.bruggeman_tortuosity.Bruggeman object at 0x13e80cd10>\n",
+      "through-cell convection <pybamm.models.submodels.convection.through_cell.no_convection.NoConvection object at 0x13e80ce50>\n",
+      "transverse convection <pybamm.models.submodels.convection.transverse.no_convection.NoConvection object at 0x13e80cb50>\n",
+      "negative interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x13e80c750>\n",
+      "positive interface <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.InverseButlerVolmer object at 0x13e81d050>\n",
+      "negative interface current <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x13e81d150>\n",
+      "positive interface current <pybamm.models.submodels.interface.inverse_kinetics.inverse_butler_volmer.CurrentForInverseButlerVolmer object at 0x13e81d250>\n",
+      "negative oxygen interface <pybamm.models.submodels.interface.kinetics.no_reaction.NoReaction object at 0x13e81d350>\n",
+      "positive oxygen interface <pybamm.models.submodels.interface.kinetics.no_reaction.NoReaction object at 0x13e81d450>\n",
+      "negative particle <pybamm.models.submodels.particle.fast_single_particle.FastSingleParticle object at 0x13dc80d10>\n",
+      "positive particle <pybamm.models.submodels.particle.fickian_single_particle.FickianSingleParticle object at 0x13e81d650>\n",
+      "negative electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x13e81d750>\n",
+      "leading-order electrolyte conductivity <pybamm.models.submodels.electrolyte_conductivity.leading_order_conductivity.LeadingOrder object at 0x13e81d850>\n",
+      "electrolyte diffusion <pybamm.models.submodels.electrolyte_diffusion.constant_concentration.ConstantConcentration object at 0x13e81d950>\n",
+      "positive electrode <pybamm.models.submodels.electrode.ohm.leading_ohm.LeadingOrder object at 0x13e81da50>\n",
+      "thermal <pybamm.models.submodels.thermal.isothermal.Isothermal object at 0x13e81db50>\n",
+      "current collector <pybamm.models.submodels.current_collector.homogeneous_current_collector.Uniform object at 0x13e81dc50>\n",
+      "negative sei <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x13e81dd50>\n",
+      "positive sei <pybamm.models.submodels.interface.sei.no_sei.NoSEI object at 0x13e81de50>\n"
      ]
     }
    ],
@@ -183,9 +178,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 7,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{}"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "model.rhs"
    ]
@@ -199,7 +205,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -215,15 +221,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "{Variable(0xed7452bed2f1969, Discharge capacity [A.h], children=[], domain=[], auxiliary_domains={}): Division(0x381d630aa1528443, /, children=['Current function [A] * 96485.33212 * Maximum concentration in negative electrode [mol.m-3] * Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m] / function (absolute)', '3600.0'], domain=[], auxiliary_domains={}),\n",
-       " Variable(0x5f6d314a2ae28139, X-averaged negative particle surface concentration, children=[], domain=['current collector'], auxiliary_domains={}): Division(-0x51e6c8f5bc31c548, /, children=['-3.0 * broadcast(Current function [A] / Typical current [A] * function (sign)) / Negative electrode thickness [m] / Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m]', 'Negative electrode surface area to volume ratio [m-1] * Negative particle radius [m]'], domain=['current collector'], auxiliary_domains={}),\n",
-       " Variable(0x25e9b81e826ecc78, X-averaged positive particle concentration, children=[], domain=['positive particle'], auxiliary_domains={'secondary': \"['current collector']\"}): Multiplication(0x2cd744f1f248da68, *, children=['-1.0 / Positive particle radius [m] ** 2.0 / Positive electrode diffusivity [m2.s-1] / 96485.33212 * Maximum concentration in negative electrode [mol.m-3] * Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m] / function (absolute)', 'div(-Positive electrode diffusivity [m2.s-1] / Positive electrode diffusivity [m2.s-1] * grad(X-averaged positive particle concentration))'], domain=['positive particle'], auxiliary_domains={'secondary': \"['current collector']\"})}"
+       "{Variable(-0x4f68261d0a00d98b, Discharge capacity [A.h], children=[], domain=[], auxiliary_domains={}): Division(-0x1fe7fec2f8387d95, /, children=['Current function [A] * 96485.33212 * Maximum concentration in negative electrode [mol.m-3] * Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m] / function (absolute)', '3600.0'], domain=[], auxiliary_domains={}),\n",
+       " Variable(0x3a35670ab4562c18, X-averaged negative particle surface concentration, children=[], domain=['current collector'], auxiliary_domains={}): Division(0x3ec3c463e3cd4a1b, /, children=[\"-3.0 * integral dx_n ['negative electrode'](broadcast(broadcast(Current function [A] / Typical current [A] * function (sign)) / Negative electrode thickness [m] / Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m]) - broadcast(0.0) + broadcast(0.0)) / Negative electrode thickness [m] / Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m]\", 'Negative electrode surface area to volume ratio [m-1] * Negative particle radius [m]'], domain=['current collector'], auxiliary_domains={}),\n",
+       " Variable(-0xa0b4a781ec8b92e, X-averaged positive particle concentration, children=[], domain=['positive particle'], auxiliary_domains={'secondary': \"['current collector']\"}): Multiplication(-0x2f5dd5cda9b83ae0, *, children=['-1.0 / Positive particle radius [m] ** 2.0 / Positive electrode diffusivity [m2.s-1] / 96485.33212 * Maximum concentration in negative electrode [mol.m-3] * Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m] / function (absolute)', 'div(-Positive electrode diffusivity [m2.s-1] / Positive electrode diffusivity [m2.s-1] * grad(X-averaged positive particle concentration))'], domain=['positive particle'], auxiliary_domains={'secondary': \"['current collector']\"})}"
       ]
      },
      "execution_count": 9,
@@ -244,7 +250,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
@@ -302,7 +308,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -320,7 +326,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -390,7 +396,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the Single Particle Model, the overpotential can be obtianed by inverting the Butler-Volmer relation, so we choose the `InverseButlerVolmer` submodel for the interface, with the \"main\" lithium-ion reaction"
+    "In the Single Particle Model, the overpotential can be obtianed by inverting the Butler-Volmer relation, so we choose the `InverseButlerVolmer` submodel for the interface, with the \"main\" lithium-ion reaction. Because of how the current is implemented, we also need to separately specify the `CurrentForInverseButlerVolmer` submodel"
    ]
   },
   {
@@ -404,14 +410,24 @@
     "] = pybamm.interface.InverseButlerVolmer(model.param, \"Negative\", \"lithium-ion main\")\n",
     "model.submodels[\n",
     "    \"positive interface\"\n",
-    "] = pybamm.interface.InverseButlerVolmer(model.param, \"Positive\", \"lithium-ion main\")\n"
+    "] = pybamm.interface.InverseButlerVolmer(model.param, \"Positive\", \"lithium-ion main\")\n",
+    "model.submodels[\n",
+    "    \"negative interface current\"\n",
+    "] = pybamm.interface.CurrentForInverseButlerVolmer(\n",
+    "    model.param, \"Negative\", \"lithium-ion main\"\n",
+    ")\n",
+    "model.submodels[\n",
+    "    \"positive interface current\"\n",
+    "] = pybamm.interface.CurrentForInverseButlerVolmer(\n",
+    "    model.param, \"Positive\", \"lithium-ion main\"\n",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Finally, for the electrolyte we assume that diffusion is infinitely fast so that the concentration is uniform, and also use the leading-order model for charge conservation, which leads to a linear variation in ionic current in the electrodes"
+    "We don't want any SEI formation in this model"
    ]
   },
   {
@@ -419,6 +435,23 @@
    "execution_count": 17,
    "metadata": {},
    "outputs": [],
+   "source": [
+    "model.submodels[\"negative sei\"] = pybamm.sei.NoSEI(model.param, \"Negative\")\n",
+    "model.submodels[\"positive sei\"] = pybamm.sei.NoSEI(model.param, \"Positive\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally, for the electrolyte we assume that diffusion is infinitely fast so that the concentration is uniform, and also use the leading-order model for charge conservation, which leads to a linear variation in ionic current in the electrodes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "model.submodels[\"electrolyte diffusion\"] = pybamm.electrolyte_diffusion.ConstantConcentration(\n",
     "    model.param\n",
@@ -437,7 +470,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -453,7 +486,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -469,7 +502,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [
     {
diff --git a/examples/scripts/sei_growth.py b/examples/scripts/calendar_ageing.py
similarity index 90%
rename from examples/scripts/sei_growth.py
rename to examples/scripts/calendar_ageing.py
index 0391bef490..bd591742aa 100644
--- a/examples/scripts/sei_growth.py
+++ b/examples/scripts/calendar_ageing.py
@@ -6,12 +6,7 @@
 models = [
     pb.lithium_ion.SPM({"sei": "reaction limited"}),
     pb.lithium_ion.SPM(
-        {
-            "sei": "reaction limited",
-            # "sei film resistance": "average",
-            "surface form": "algebraic",
-        },
-        name="Algebraic SPM",
+        {"sei": "reaction limited", "surface form": "algebraic"}, name="Algebraic SPM",
     ),
     pb.lithium_ion.DFN({"sei": "reaction limited"}),
 ]
diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py
index 9c3258f1fd..b81aaf8b5a 100644
--- a/examples/scripts/compare_lead_acid.py
+++ b/examples/scripts/compare_lead_acid.py
@@ -54,6 +54,10 @@
     "Porosity",
     "Electrolyte potential [V]",
     "Terminal voltage [V]",
+    "Negative electrode reaction overpotential",
+    "Positive electrode reaction overpotential",
+    "Sum of interfacial current densities",
+    "Sum of electrolyte reaction source terms",
 ]
 plot = pybamm.QuickPlot(solutions, output_variables, linestyles=[":", "--", "-"])
 plot.dynamic_plot()
diff --git a/examples/scripts/custom_model.py b/examples/scripts/custom_model.py
index 8d0879ae0e..ae8c042723 100644
--- a/examples/scripts/custom_model.py
+++ b/examples/scripts/custom_model.py
@@ -34,12 +34,24 @@
 model.submodels["positive interface"] = pybamm.interface.InverseButlerVolmer(
     model.param, "Positive", "lithium-ion main"
 )
+model.submodels[
+    "negative interface current"
+] = pybamm.interface.CurrentForInverseButlerVolmer(
+    model.param, "Negative", "lithium-ion main"
+)
+model.submodels[
+    "positive interface current"
+] = pybamm.interface.CurrentForInverseButlerVolmer(
+    model.param, "Positive", "lithium-ion main"
+)
 model.submodels[
     "electrolyte diffusion"
 ] = pybamm.electrolyte_diffusion.ConstantConcentration(model.param)
 model.submodels[
     "electrolyte conductivity"
 ] = pybamm.electrolyte_conductivity.LeadingOrder(model.param)
+model.submodels["negative sei"] = pybamm.sei.NoSEI(model.param, "Negative")
+model.submodels["positive sei"] = pybamm.sei.NoSEI(model.param, "Positive")
 
 # build model
 model.build_model()
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index 97cafcb6e3..0bf88228a3 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -144,12 +144,12 @@ def set_interfacial_submodel(self):
             self.submodels[
                 "negative interface current"
             ] = pybamm.interface.CurrentForInverseButlerVolmer(
-                self.param, "Negative", "lithium-ion main"
+                self.param, "Negative", "lead-acid main"
             )
             self.submodels[
                 "positive interface current"
             ] = pybamm.interface.CurrentForInverseButlerVolmer(
-                self.param, "Positive", "lithium-ion main"
+                self.param, "Positive", "lead-acid main"
             )
         else:
             self.submodels[
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 fe084c068d..8b1abcdb43 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -23,7 +23,7 @@ class InverseButlerVolmer(BaseInterface):
         A dictionary of options to be passed to the model. In this case "sei film
         resistance" is the important option. See :class:`pybamm.BaseBatteryModel`
 
-    **Extends:** :class:`pybamm.interface.kinetics.ButlerVolmer`
+    **Extends:** :class:`pybamm.interface.BaseInterface`
 
     """
 
@@ -118,7 +118,7 @@ class CurrentForInverseButlerVolmer(BaseInterface):
     reaction : str
         The name of the reaction being implemented
 
-    **Extends:** :class:`pybamm.interface.kinetics.ButlerVolmer`
+    **Extends:** :class:`pybamm.interface.BaseInterface`
 
     """
 
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index e1f44d97e9..5b09b36d46 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -530,6 +530,14 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.j_p_av = solution[
             "X-averaged positive electrode interfacial current density"
         ]
+        self.j_n_sei = solution["Negative electrode sei interfacial current density"]
+        self.j_p_sei = solution["Positive electrode sei interfacial current density"]
+        self.j_n_sei_av = solution[
+            "X-averaged negative electrode sei interfacial current density"
+        ]
+        self.j_p_sei_av = solution[
+            "X-averaged positive electrode sei interfacial current density"
+        ]
 
         self.j0_n = solution["Negative electrode exchange current density"]
         self.j0_p = solution["Positive electrode exchange current density"]
@@ -543,10 +551,14 @@ def test_interfacial_current_average(self):
         """Test that average of the interfacial current density is equal to the true
         value."""
         np.testing.assert_array_almost_equal(
-            self.j_n_av(self.t), self.i_cell / self.l_n, decimal=4
+            self.j_n_av(self.t) + self.j_n_sei_av(self.t),
+            self.i_cell / self.l_n,
+            decimal=4,
         )
         np.testing.assert_array_almost_equal(
-            self.j_p_av(self.t), -self.i_cell / self.l_p, decimal=4
+            self.j_p_av(self.t) + self.j_p_sei_av(self.t),
+            -self.i_cell / self.l_p,
+            decimal=4,
         )
 
     def test_conservation(self):

From c6888164d6457d9f0f3388d26e8cdc49e70f036a Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 14 May 2020 09:50:49 -0400
Subject: [PATCH 06/15] #984 coverage

---
 .../test_full_battery_models/test_lithium_ion/test_dfn.py   | 6 ++++++
 .../test_full_battery_models/test_lithium_ion/test_dfn.py   | 5 +++++
 .../test_interface/test_kinetics/test_tafel.py              | 4 +++-
 3 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 8bfe950786..313b338daa 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -124,6 +124,12 @@ def positive_distribution(x):
 
 
 class TestDFNWithSEI(unittest.TestCase):
+    def test_well_posed_constant(self):
+        options = {"sei": "constant"}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
     def test_well_posed_reaction_limited(self):
         options = {"sei": "reaction limited"}
         model = pybamm.lithium_ion.DFN(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index a4d201b8a8..53c7919b3a 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -119,6 +119,11 @@ def test_surface_form_algebraic(self):
 
 
 class TestDFNWithSEI(unittest.TestCase):
+    def test_well_posed_constant(self):
+        options = {"sei": "constant"}
+        model = pybamm.lithium_ion.DFN(options)
+        model.check_well_posedness()
+
     def test_well_posed_reaction_limited(self):
         options = {"sei": "reaction limited"}
         model = pybamm.lithium_ion.DFN(options)
diff --git a/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_tafel.py b/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_tafel.py
index bcbec88e76..2eacdf6be5 100644
--- a/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_tafel.py
+++ b/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_tafel.py
@@ -8,6 +8,7 @@
 
 class TestTafel(unittest.TestCase):
     def test_public_function(self):
+        # Test forward Tafel on the negative
         param = pybamm.standard_parameters_lithium_ion
 
         a_n = pybamm.FullBroadcast(
@@ -31,6 +32,7 @@ def test_public_function(self):
 
         std_tests.test_all()
 
+        # Test backward Tafel on the positive
         variables = {
             "Current collector current density": a,
             "Positive electrode potential": a_p,
@@ -56,7 +58,7 @@ def test_public_function(self):
             "Sum of x-averaged negative electrode interfacial current densities": 0,
             "Sum of x-averaged positive electrode interfacial current densities": 0,
         }
-        submodel = pybamm.interface.ForwardTafel(param, "Positive", "lithium-ion main")
+        submodel = pybamm.interface.BackwardTafel(param, "Positive", "lithium-ion main")
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 

From 5a55dbdcab24f81e4b1be5054d54a9fafdf64b63 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 14 May 2020 13:42:14 -0400
Subject: [PATCH 07/15] #984 fix scalings

---
 examples/scripts/calendar_ageing.py              |  8 ++++++--
 examples/scripts/create-model.py                 |  4 ++--
 .../lithium-ion/seis/Example/parameters.csv      |  2 +-
 .../models/submodels/interface/base_interface.py |  8 ++++----
 .../inverse_kinetics/inverse_butler_volmer.py    |  9 ++++-----
 .../models/submodels/interface/sei/base_sei.py   | 16 ++++++++++------
 .../submodels/interface/sei/reaction_limited.py  |  4 +++-
 7 files changed, 30 insertions(+), 21 deletions(-)

diff --git a/examples/scripts/calendar_ageing.py b/examples/scripts/calendar_ageing.py
index bd591742aa..b0bdf731b7 100644
--- a/examples/scripts/calendar_ageing.py
+++ b/examples/scripts/calendar_ageing.py
@@ -21,7 +21,7 @@
 
     solver = pb.CasadiSolver(mode="fast")
 
-    years = 3
+    years = 30
     days = years * 365
     hours = days * 24
     minutes = hours * 60
@@ -42,7 +42,7 @@
         "X-averaged total negative electrode sei thickness [m]",
         "X-averaged total negative electrode sei thickness",
         "X-averaged negative electrode sei concentration [mol.m-3]",
-        "Loss of lithium to negative electrode sei [mols]",
+        "Loss of lithium to negative electrode sei [mol]",
         [
             "Negative electrode sei interfacial current density [A.m-2]",
             "Negative electrode interfacial current density [A.m-2]",
@@ -51,6 +51,10 @@
             "X-averaged negative electrode sei interfacial current density [A.m-2]",
             "X-averaged negative electrode interfacial current density [A.m-2]",
         ],
+        [
+            "X-averaged negative electrode sei interfacial current density",
+            "X-averaged negative electrode interfacial current density",
+        ],
         "Sum of x-averaged negative electrode interfacial current densities",
         "Sum of negative electrode interfacial current densities",
         "X-averaged electrolyte concentration",
diff --git a/examples/scripts/create-model.py b/examples/scripts/create-model.py
index 73b9161351..3b1d8c3467 100644
--- a/examples/scripts/create-model.py
+++ b/examples/scripts/create-model.py
@@ -67,8 +67,8 @@ def D(cc):
     "SEI growth rate": dLdt,
     "Solvent concentration": c,
     "SEI thickness [m]": L_0_dim * L,
-    "SEI growth rate [m/s]": (D_dim(c_inf_dim) / L_0_dim) * dLdt,
-    "Solvent concentration [mols/m^3]": c_inf_dim * c,
+    "SEI growth rate [m.s-1]": (D_dim(c_inf_dim) / L_0_dim) * dLdt,
+    "Solvent concentration [mol.m-3]": c_inf_dim * c,
 }
 
 "--------------------------------------------------------------------------------------"
diff --git a/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv b/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
index 638b8d8aa3..63fe096c5f 100644
--- a/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
+++ b/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
@@ -6,7 +6,7 @@ Inner SEI reaction proportion,0.5,,
 Inner SEI partial molar volume [m3.mol-1],3e-6, Guess,
 Outer SEI partial molar volume [m3.mol-1],3e-6, Guess,
 SEI reaction exchange current density [A.m-2],1.5E-7, Guess,
-SEI resistivity [Ohm.m],1, Ramadass,
+SEI resistivity [Ohm.m],1e6,Guess,check Ramadass
 Outer SEI solvent diffusivity [m2.s-1],2.5E-22, Single paper,
 Bulk solvent concentration [mol.m-3],2.636E3, Ploehn paper,
 Ratio of inner and outer SEI exchange current densities,1, Assume same,
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index 75f9346cdf..408d0fbfbd 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -35,9 +35,9 @@ def __init__(self, param, domain, reaction):
         elif reaction == "lithium-ion oxygen":
             self.reaction_name = " oxygen"
             self.Reaction_icd = "Oxygen interfacial current density"
-        elif reaction == "sei":
-            self.reaction_name = " sei"
-            self.Reaction_icd = "Sei interfacial current density"
+        elif reaction == "scaled sei":
+            self.reaction_name = " scaled sei"
+            self.Reaction_icd = "Scaled sei interfacial current density"
         self.reaction = reaction
 
     def _get_exchange_current_density(self, variables):
@@ -175,7 +175,7 @@ def _get_number_of_electrons_in_reaction(self):
 
     def _get_electrolyte_reaction_signed_stoichiometry(self):
         "Returns the number of electrons in the reaction"
-        if self.reaction in ["lithium-ion main", "sei"]:
+        if self.reaction in ["lithium-ion main", "scaled sei"]:
             # Both the main reaction current contribute to the electrolyte reaction
             # current
             return pybamm.Scalar(1), pybamm.Scalar(1)
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 8b1abcdb43..6c78e32891 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -46,10 +46,7 @@ def get_coupled_variables(self, variables):
                 j_tot_av, [self.domain.lower() + " electrode"]
             )
 
-        if self.domain == "Negative":
-            ne = self.param.ne_n
-        elif self.domain == "Positive":
-            ne = self.param.ne_p
+        ne = self._get_number_of_electrons_in_reaction()
         # Note: T must have the same domain as j0 and eta_r
         if j0.domain in ["current collector", ["current collector"]]:
             T = variables["X-averaged cell temperature"]
@@ -131,7 +128,9 @@ def get_coupled_variables(self, variables):
             + self.domain.lower()
             + " electrode total interfacial current density"
         ]
-        j_sei = variables[self.domain + " electrode sei interfacial current density"]
+        j_sei = variables[
+            self.domain + " electrode scaled sei interfacial current density"
+        ]
         j = j_tot - j_sei
 
         variables.update(self._get_standard_interfacial_current_variables(j))
diff --git a/pybamm/models/submodels/interface/sei/base_sei.py b/pybamm/models/submodels/interface/sei/base_sei.py
index aee9408d06..2feba30c61 100644
--- a/pybamm/models/submodels/interface/sei/base_sei.py
+++ b/pybamm/models/submodels/interface/sei/base_sei.py
@@ -19,7 +19,7 @@ class BaseModel(BaseInterface):
     """
 
     def __init__(self, param, domain):
-        reaction = "sei"
+        reaction = "scaled sei"
         super().__init__(param, domain, reaction)
 
     def _get_standard_thickness_variables(self, L_inner, L_outer):
@@ -95,7 +95,7 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
             + " sei concentration [mol.m-3]": n_outer_av * n_outer_scale,
             self.domain + " sei concentration [mol.m-3]": n_SEI * n_scale,
             "X-averaged " + domain + " sei concentration [mol.m-3]": n_SEI_av * n_scale,
-            "Loss of lithium to " + domain + " sei [mols]": Q_sei * n_scale,
+            "Loss of lithium to " + domain + " sei [mol]": Q_sei * n_scale,
         }
 
         return variables
@@ -121,13 +121,16 @@ def _get_standard_reaction_variables(self, j_inner, j_outer):
         # by parameter values in general
         if isinstance(self, pybamm.sei.NoSEI):
             j_scale = 1
-            Gamma_SEI_n = 1
+            Gamma_SEI = 1
         else:
             sp = pybamm.sei_parameters
             j_scale = (
                 sp.F * sp.L_sei_0_dim / sp.V_bar_inner_dimensional / sp.tau_discharge
             )
-            Gamma_SEI_n = sp.Gamma_SEI_n
+            if self.domain == "Negative":
+                Gamma_SEI = sp.Gamma_SEI_n
+            elif self.domain == "Positive":
+                Gamma_SEI = sp.Gamma_SEI_p
         j_i_av = pybamm.x_average(j_inner)
         j_o_av = pybamm.x_average(j_outer)
 
@@ -160,9 +163,10 @@ def _get_standard_reaction_variables(self, j_inner, j_outer):
             "X-averaged "
             + domain
             + " sei interfacial current density [A.m-2]": j_sei_av * j_scale,
-            "Scaled "
+            Domain + " scaled sei interfacial current density": j_sei * Gamma_SEI,
+            "X-averaged "
             + domain
-            + " sei interfacial current density": j_sei * Gamma_SEI_n,
+            + " scaled sei interfacial current density": j_sei_av * Gamma_SEI,
         }
 
         return variables
diff --git a/pybamm/models/submodels/interface/sei/reaction_limited.py b/pybamm/models/submodels/interface/sei/reaction_limited.py
index 89a250949b..69aa589e8c 100644
--- a/pybamm/models/submodels/interface/sei/reaction_limited.py
+++ b/pybamm/models/submodels/interface/sei/reaction_limited.py
@@ -54,7 +54,9 @@ def get_coupled_variables(self, variables):
         # alpha = pybamm.sei_parameters.alpha
 
         # need to revise for thermal case
-        j_sei = -(1 / C_sei) * pybamm.exp(-(phi_s_n - phi_e_n - j * L_sei * R_sei))
+        j_sei = -(1 / C_sei) * pybamm.exp(
+            -0.5 * (phi_s_n - phi_e_n - j * L_sei * R_sei)
+        )
 
         j_inner = alpha * j_sei
         j_outer = (1 - alpha) * j_sei

From daf49edde3a18898ceac83eba8e0501987b0b93c Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 14 May 2020 13:49:31 -0400
Subject: [PATCH 08/15] #984 reformat interface docs

---
 .../first_order_kinetics.rst                  |  0
 .../interface/first_order_kinetics/index.rst  |  8 -------
 .../inverse_first_order_kinetics.rst          |  5 ----
 .../models/submodels/interface/index.rst      |  8 +++----
 ...butler_volmer.rst => inverse_kinetics.rst} |  4 ++--
 .../interface/inverse_kinetics/index.rst      |  7 ------
 .../models/submodels/interface/kinetics.rst   | 17 ++++++++++++++
 .../interface/kinetics/base_kinetics.rst      |  5 ----
 .../interface/kinetics/butler_volmer.rst      |  5 ----
 .../submodels/interface/kinetics/index.rst    | 10 --------
 .../interface/kinetics/no_reaction.rst        |  5 ----
 .../submodels/interface/kinetics/tafel.rst    |  8 -------
 .../source/models/submodels/interface/sei.rst | 23 +++++++++++++++++++
 .../submodels/interface/sei/base_sei.rst      |  6 -----
 .../sei/electron_migration_limited.rst        |  6 -----
 .../models/submodels/interface/sei/index.rst  | 12 ----------
 .../sei/interstitial_diffusion_limited.rst    |  6 -----
 .../models/submodels/interface/sei/no_sei.rst |  6 -----
 .../interface/sei/reaction_limited.rst        |  6 -----
 .../sei/solvent_diffusion_limited.rst         |  6 -----
 20 files changed, 46 insertions(+), 107 deletions(-)
 rename docs/source/models/submodels/interface/{first_order_kinetics => }/first_order_kinetics.rst (100%)
 delete mode 100644 docs/source/models/submodels/interface/first_order_kinetics/index.rst
 delete mode 100644 docs/source/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.rst
 rename docs/source/models/submodels/interface/{inverse_kinetics/inverse_butler_volmer.rst => inverse_kinetics.rst} (65%)
 delete mode 100644 docs/source/models/submodels/interface/inverse_kinetics/index.rst
 create mode 100644 docs/source/models/submodels/interface/kinetics.rst
 delete mode 100644 docs/source/models/submodels/interface/kinetics/base_kinetics.rst
 delete mode 100644 docs/source/models/submodels/interface/kinetics/butler_volmer.rst
 delete mode 100644 docs/source/models/submodels/interface/kinetics/index.rst
 delete mode 100644 docs/source/models/submodels/interface/kinetics/no_reaction.rst
 delete mode 100644 docs/source/models/submodels/interface/kinetics/tafel.rst
 create mode 100644 docs/source/models/submodels/interface/sei.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/base_sei.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/electron_migration_limited.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/index.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/interstitial_diffusion_limited.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/no_sei.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/reaction_limited.rst
 delete mode 100644 docs/source/models/submodels/interface/sei/solvent_diffusion_limited.rst

diff --git a/docs/source/models/submodels/interface/first_order_kinetics/first_order_kinetics.rst b/docs/source/models/submodels/interface/first_order_kinetics.rst
similarity index 100%
rename from docs/source/models/submodels/interface/first_order_kinetics/first_order_kinetics.rst
rename to docs/source/models/submodels/interface/first_order_kinetics.rst
diff --git a/docs/source/models/submodels/interface/first_order_kinetics/index.rst b/docs/source/models/submodels/interface/first_order_kinetics/index.rst
deleted file mode 100644
index fab79b9b44..0000000000
--- a/docs/source/models/submodels/interface/first_order_kinetics/index.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-First-order Kinetics
-====================
-
-.. toctree::
-  :maxdepth: 1
-
-  inverse_first_order_kinetics
-  first_order_kinetics
diff --git a/docs/source/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.rst b/docs/source/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.rst
deleted file mode 100644
index 58e1fb83b5..0000000000
--- a/docs/source/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Base Inverse First-order Kinetics
-=================================
-
-.. autoclass:: pybamm.interface.InverseFirstOrderKinetics
-    :members:
diff --git a/docs/source/models/submodels/interface/index.rst b/docs/source/models/submodels/interface/index.rst
index c2fb250218..f616791487 100644
--- a/docs/source/models/submodels/interface/index.rst
+++ b/docs/source/models/submodels/interface/index.rst
@@ -5,8 +5,8 @@ Interface
   :maxdepth: 1
 
   base_interface
-  kinetics/index
-  inverse_kinetics/index
-  first_order_kinetics/index
+  kinetics
+  inverse_kinetics
+  first_order_kinetics
   diffusion_limited
-  sei/index
+  sei
diff --git a/docs/source/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.rst b/docs/source/models/submodels/interface/inverse_kinetics.rst
similarity index 65%
rename from docs/source/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.rst
rename to docs/source/models/submodels/interface/inverse_kinetics.rst
index ef1053594e..5318109f29 100644
--- a/docs/source/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.rst
+++ b/docs/source/models/submodels/interface/inverse_kinetics.rst
@@ -1,5 +1,5 @@
-Inverse Butler-Volmer
-=====================
+Inverse Kinetics
+================
 
 .. autoclass:: pybamm.interface.inverse_kinetics.InverseButlerVolmer
     :members:
diff --git a/docs/source/models/submodels/interface/inverse_kinetics/index.rst b/docs/source/models/submodels/interface/inverse_kinetics/index.rst
deleted file mode 100644
index d36fa61109..0000000000
--- a/docs/source/models/submodels/interface/inverse_kinetics/index.rst
+++ /dev/null
@@ -1,7 +0,0 @@
-Inverse Interface Kinetics
-==========================
-
-.. toctree::
-  :maxdepth: 1
-
-  inverse_butler_volmer
diff --git a/docs/source/models/submodels/interface/kinetics.rst b/docs/source/models/submodels/interface/kinetics.rst
new file mode 100644
index 0000000000..7ed998506e
--- /dev/null
+++ b/docs/source/models/submodels/interface/kinetics.rst
@@ -0,0 +1,17 @@
+Kinetics
+========
+
+.. autoclass:: pybamm.interface.BaseKinetics
+    :members:
+
+.. autoclass:: pybamm.interface.ButlerVolmer
+    :members:
+
+.. autoclass:: pybamm.interface.NoReaction
+    :members:
+
+.. autoclass:: pybamm.interface.ForwardTafel
+    :members:
+
+.. autoclass:: pybamm.interface.BackwardTafel
+    :members:
diff --git a/docs/source/models/submodels/interface/kinetics/base_kinetics.rst b/docs/source/models/submodels/interface/kinetics/base_kinetics.rst
deleted file mode 100644
index b91c17f810..0000000000
--- a/docs/source/models/submodels/interface/kinetics/base_kinetics.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Base Kinetics
-=============
-
-.. autoclass:: pybamm.interface.BaseKinetics
-    :members:
diff --git a/docs/source/models/submodels/interface/kinetics/butler_volmer.rst b/docs/source/models/submodels/interface/kinetics/butler_volmer.rst
deleted file mode 100644
index 8d456d9bf8..0000000000
--- a/docs/source/models/submodels/interface/kinetics/butler_volmer.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Butler-Volmer
-=============
-
-.. autoclass:: pybamm.interface.ButlerVolmer
-    :members:
diff --git a/docs/source/models/submodels/interface/kinetics/index.rst b/docs/source/models/submodels/interface/kinetics/index.rst
deleted file mode 100644
index 3dc0fcc859..0000000000
--- a/docs/source/models/submodels/interface/kinetics/index.rst
+++ /dev/null
@@ -1,10 +0,0 @@
-Interface Kinetics
-==================
-
-.. toctree::
-  :maxdepth: 1
-
-  base_kinetics
-  butler_volmer
-  no_reaction
-  tafel
diff --git a/docs/source/models/submodels/interface/kinetics/no_reaction.rst b/docs/source/models/submodels/interface/kinetics/no_reaction.rst
deleted file mode 100644
index 3d1e91455d..0000000000
--- a/docs/source/models/submodels/interface/kinetics/no_reaction.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-No Reaction
-===========
-
-.. autoclass:: pybamm.interface.NoReaction
-    :members:
diff --git a/docs/source/models/submodels/interface/kinetics/tafel.rst b/docs/source/models/submodels/interface/kinetics/tafel.rst
deleted file mode 100644
index dcb425214b..0000000000
--- a/docs/source/models/submodels/interface/kinetics/tafel.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-Tafel
-=====
-
-.. autoclass:: pybamm.interface.ForwardTafel
-    :members:
-
-.. autoclass:: pybamm.interface.BackwardTafel
-    :members:
diff --git a/docs/source/models/submodels/interface/sei.rst b/docs/source/models/submodels/interface/sei.rst
new file mode 100644
index 0000000000..a036263f27
--- /dev/null
+++ b/docs/source/models/submodels/interface/sei.rst
@@ -0,0 +1,23 @@
+SEI models
+==========
+
+.. autoclass:: pybamm.sei.BaseModel
+    :members:
+
+.. autoclass:: pybamm.sei.ConstantSEI
+    :members:
+
+.. autoclass:: pybamm.sei.ElectronMigrationLimited
+    :members:
+
+.. autoclass:: pybamm.sei.InterstitialDiffusionLimited
+    :members:
+
+.. autoclass:: pybamm.sei.NoSEI
+    :members:
+
+.. autoclass:: pybamm.sei.ReactionLimited
+    :members:
+
+.. autoclass:: pybamm.sei.SolventDiffusionLimited
+    :members:
\ No newline at end of file
diff --git a/docs/source/models/submodels/interface/sei/base_sei.rst b/docs/source/models/submodels/interface/sei/base_sei.rst
deleted file mode 100644
index 77e52fa595..0000000000
--- a/docs/source/models/submodels/interface/sei/base_sei.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Base SEI
-========
-
-.. autoclass:: pybamm.sei.BaseModel
-    :members:
-
diff --git a/docs/source/models/submodels/interface/sei/electron_migration_limited.rst b/docs/source/models/submodels/interface/sei/electron_migration_limited.rst
deleted file mode 100644
index 8d969803d3..0000000000
--- a/docs/source/models/submodels/interface/sei/electron_migration_limited.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Electron Migration Limited
-==========================
-
-.. autoclass:: pybamm.sei.ElectronMigrationLimited
-    :members:
-
diff --git a/docs/source/models/submodels/interface/sei/index.rst b/docs/source/models/submodels/interface/sei/index.rst
deleted file mode 100644
index d7a138d47d..0000000000
--- a/docs/source/models/submodels/interface/sei/index.rst
+++ /dev/null
@@ -1,12 +0,0 @@
-SEI
-===
-
-.. toctree::
-  :maxdepth: 1
-
-  base_sei
-  electron_migration_limited
-  interstitial_diffusion_limited
-  no_sei
-  reaction_limited
-  solvent_diffusion_limited
diff --git a/docs/source/models/submodels/interface/sei/interstitial_diffusion_limited.rst b/docs/source/models/submodels/interface/sei/interstitial_diffusion_limited.rst
deleted file mode 100644
index 244b2dbfb8..0000000000
--- a/docs/source/models/submodels/interface/sei/interstitial_diffusion_limited.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Interstitial Diffusion Limited
-==============================
-
-.. autoclass:: pybamm.sei.InterstitialDiffusionLimited
-    :members:
-
diff --git a/docs/source/models/submodels/interface/sei/no_sei.rst b/docs/source/models/submodels/interface/sei/no_sei.rst
deleted file mode 100644
index 3bda8291e0..0000000000
--- a/docs/source/models/submodels/interface/sei/no_sei.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-No SEI
-======
-
-.. autoclass:: pybamm.sei.NoSEI
-    :members:
-
diff --git a/docs/source/models/submodels/interface/sei/reaction_limited.rst b/docs/source/models/submodels/interface/sei/reaction_limited.rst
deleted file mode 100644
index 09ee42175d..0000000000
--- a/docs/source/models/submodels/interface/sei/reaction_limited.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Reaction Limited
-================
-
-.. autoclass:: pybamm.sei.ReactionLimited
-    :members:
-
diff --git a/docs/source/models/submodels/interface/sei/solvent_diffusion_limited.rst b/docs/source/models/submodels/interface/sei/solvent_diffusion_limited.rst
deleted file mode 100644
index 07ec74e71d..0000000000
--- a/docs/source/models/submodels/interface/sei/solvent_diffusion_limited.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Solvent Diffusion Limited
-=========================
-
-.. autoclass:: pybamm.sei.SolventDiffusionLimited
-    :members:
-

From e574b4f6ef05f79aecb8dbaa83548a84b60ce83c Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 14 May 2020 13:50:20 -0400
Subject: [PATCH 09/15] #984 add unsaved file

---
 .../source/models/submodels/interface/first_order_kinetics.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/docs/source/models/submodels/interface/first_order_kinetics.rst b/docs/source/models/submodels/interface/first_order_kinetics.rst
index 32522d6205..874b1a9a87 100644
--- a/docs/source/models/submodels/interface/first_order_kinetics.rst
+++ b/docs/source/models/submodels/interface/first_order_kinetics.rst
@@ -3,3 +3,6 @@ First-order Kinetics
 
 .. autoclass:: pybamm.interface.FirstOrderKinetics
     :members:
+
+.. autoclass:: pybamm.interface.InverseFirstOrderKinetics
+    :members:

From aaa61dd0eb31dfc21172e1c97a5b1f8a45b5ee3c Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 14 May 2020 14:59:37 -0400
Subject: [PATCH 10/15] #984 fix models for output tests

---
 .../full_battery_models/base_battery_model.py |  18 ++
 .../leading_order_conductivity.py             |   8 +-
 .../submodels/interface/base_interface.py     |  24 ++
 .../submodels/interface/diffusion_limited.py  |   4 +
 .../first_order_kinetics.py                   |   4 +
 .../inverse_first_order_kinetics.py           |   5 +
 .../inverse_kinetics/inverse_butler_volmer.py |   7 +-
 .../interface/kinetics/base_kinetics.py       |  12 +-
 .../test_models/standard_output_tests.py      |  17 +-
 .../test_lithium_ion/test_spm.py              | 214 +++++++++---------
 10 files changed, 194 insertions(+), 119 deletions(-)

diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index cba817bba0..d8bebe237f 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -658,6 +658,22 @@ def set_voltage_variables(self):
         eta_r_av = eta_r_p_av - eta_r_n_av
         eta_r_av_dim = eta_r_p_av_dim - eta_r_n_av_dim
 
+        # SEI film overpotential
+        eta_sei_n_av = self.variables[
+            "X-averaged negative electrode sei film overpotential"
+        ]
+        eta_sei_p_av = self.variables[
+            "X-averaged positive electrode sei film overpotential"
+        ]
+        eta_sei_n_av_dim = self.variables[
+            "X-averaged negative electrode sei film overpotential [V]"
+        ]
+        eta_sei_p_av_dim = self.variables[
+            "X-averaged positive electrode sei film overpotential [V]"
+        ]
+        eta_sei_av = eta_sei_n_av + eta_sei_p_av
+        eta_sei_av_dim = eta_sei_n_av_dim + eta_sei_p_av_dim
+
         # TODO: add current collector losses to the voltage in 3D
 
         self.variables.update(
@@ -668,6 +684,8 @@ def set_voltage_variables(self):
                 "Measured open circuit voltage [V]": ocv_dim,
                 "X-averaged reaction overpotential": eta_r_av,
                 "X-averaged reaction overpotential [V]": eta_r_av_dim,
+                "X-averaged sei film overpotential": eta_sei_av,
+                "X-averaged sei film overpotential [V]": eta_sei_av_dim,
                 "X-averaged solid phase ohmic losses": delta_phi_s_av,
                 "X-averaged solid phase ohmic losses [V]": delta_phi_s_av_dim,
             }
diff --git a/pybamm/models/submodels/electrolyte_conductivity/leading_order_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/leading_order_conductivity.py
index 17962268af..437c94e4ee 100644
--- a/pybamm/models/submodels/electrolyte_conductivity/leading_order_conductivity.py
+++ b/pybamm/models/submodels/electrolyte_conductivity/leading_order_conductivity.py
@@ -26,10 +26,12 @@ def __init__(self, param, domain=None):
         super().__init__(param, domain)
 
     def get_coupled_variables(self, variables):
-        ocp_n_av = variables["X-averaged negative electrode open circuit potential"]
-        eta_r_n_av = variables["X-averaged negative electrode reaction overpotential"]
+        # delta_phi = phi_s - phi_e
+        delta_phi_n_av = variables[
+            "X-averaged negative electrode surface potential difference"
+        ]
         phi_s_n_av = variables["X-averaged negative electrode potential"]
-        phi_e_av = phi_s_n_av - eta_r_n_av - ocp_n_av
+        phi_e_av = phi_s_n_av - delta_phi_n_av
         return self._get_coupled_variables_from_potential(variables, phi_e_av)
 
     def _get_coupled_variables_from_potential(self, variables, phi_e_av):
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index 408d0fbfbd..cb84e63b98 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -481,6 +481,30 @@ def _get_standard_overpotential_variables(self, eta_r):
 
         return variables
 
+    def _get_standard_sei_film_overpotential_variables(self, eta_sei):
+
+        pot_scale = self.param.potential_scale
+        # Average, and broadcast if necessary
+        eta_sei_av = pybamm.x_average(eta_sei)
+        if eta_sei.domain == []:
+            eta_sei = pybamm.FullBroadcast(
+                eta_sei, self.domain_for_broadcast, "current collector"
+            )
+        elif eta_sei.domain == ["current collector"]:
+            eta_sei = pybamm.PrimaryBroadcast(eta_sei, self.domain_for_broadcast)
+
+        domain = self.domain.lower() + " electrode"
+        variables = {
+            self.domain + " electrode sei film overpotential": eta_sei,
+            "X-averaged " + domain + " sei film overpotential": eta_sei_av,
+            self.domain + " electrode sei film overpotential [V]": eta_sei * pot_scale,
+            "X-averaged "
+            + domain
+            + " sei film overpotential [V]": eta_sei_av * pot_scale,
+        }
+
+        return variables
+
     def _get_standard_surface_potential_difference_variables(self, delta_phi):
 
         if self.domain == "Negative":
diff --git a/pybamm/models/submodels/interface/diffusion_limited.py b/pybamm/models/submodels/interface/diffusion_limited.py
index 9a7f862b4a..6de23c2789 100644
--- a/pybamm/models/submodels/interface/diffusion_limited.py
+++ b/pybamm/models/submodels/interface/diffusion_limited.py
@@ -55,6 +55,10 @@ def get_coupled_variables(self, variables):
         variables.update(self._get_standard_overpotential_variables(eta_r))
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
+        # No SEI film resistance in this model
+        eta_sei = pybamm.Scalar(0)
+        variables.update(self._get_standard_sei_film_overpotential_variables(eta_sei))
+
         if (
             "Negative electrode" + self.reaction_name + " interfacial current density"
             in variables
diff --git a/pybamm/models/submodels/interface/first_order_kinetics/first_order_kinetics.py b/pybamm/models/submodels/interface/first_order_kinetics/first_order_kinetics.py
index 264fd15c0c..f134393335 100644
--- a/pybamm/models/submodels/interface/first_order_kinetics/first_order_kinetics.py
+++ b/pybamm/models/submodels/interface/first_order_kinetics/first_order_kinetics.py
@@ -67,6 +67,10 @@ def get_coupled_variables(self, variables):
         variables.update(self._get_standard_overpotential_variables(eta_r))
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
+        # SEI film resistance not implemented in this model
+        eta_sei = pybamm.Scalar(0)
+        variables.update(self._get_standard_sei_film_overpotential_variables(eta_sei))
+
         # Add first-order averages
         j_1_bar = dj_dc_0 * pybamm.x_average(c_e_1) + dj_ddeltaphi_0 * pybamm.x_average(
             delta_phi_1
diff --git a/pybamm/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.py b/pybamm/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.py
index f931b66d8e..592303e7d0 100644
--- a/pybamm/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.py
+++ b/pybamm/models/submodels/interface/first_order_kinetics/inverse_first_order_kinetics.py
@@ -1,6 +1,7 @@
 #
 # First-order Butler-Volmer kinetics
 #
+import pybamm
 from ..base_interface import BaseInterface
 
 
@@ -75,4 +76,8 @@ def get_coupled_variables(self, variables):
             self._get_standard_surface_potential_difference_variables(delta_phi)
         )
 
+        # SEI film resistance not implemented in this model
+        eta_sei = pybamm.Scalar(0)
+        variables.update(self._get_standard_sei_film_overpotential_variables(eta_sei))
+
         return variables
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 6c78e32891..6531adfe24 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -66,10 +66,12 @@ def get_coupled_variables(self, variables):
             L_sei = variables[
                 "Total " + self.domain.lower() + " electrode sei thickness"
             ]
-            delta_phi = eta_r + ocp + j_tot * L_sei * pybamm.sei_parameters.R_sei
+            eta_sei = -j_tot * L_sei * pybamm.sei_parameters.R_sei
         # Without SEI resistance
         else:
-            delta_phi = eta_r + ocp
+            eta_sei = pybamm.Scalar(0)
+
+        delta_phi = eta_r + ocp - eta_sei
 
         variables.update(
             self._get_standard_total_interfacial_current_variables(j_tot_av)
@@ -79,6 +81,7 @@ def get_coupled_variables(self, variables):
         variables.update(
             self._get_standard_surface_potential_difference_variables(delta_phi)
         )
+        variables.update(self._get_standard_sei_film_overpotential_variables(eta_sei))
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
         return variables
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index f6c79ee968..19b41ba60c 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -79,12 +79,15 @@ def get_coupled_variables(self, variables):
             j = variables[
                 self.domain + " electrode interfacial current density variable"
             ]
-            eta_r -= j * L_sei * pybamm.sei_parameters.R_sei
+            eta_sei = -j * L_sei * pybamm.sei_parameters.R_sei
         elif self.options["sei film resistance"] == "average":
             L_sei = variables[
                 "Total " + self.domain.lower() + " electrode sei thickness"
             ]
-            eta_r -= j_tot_av * L_sei * pybamm.sei_parameters.R_sei
+            eta_sei = -j_tot_av * L_sei * pybamm.sei_parameters.R_sei
+        else:
+            eta_sei = pybamm.Scalar(0)
+        eta_r += eta_sei
 
         # Get number of electrons in reaction
         ne = self._get_number_of_electrons_in_reaction()
@@ -107,6 +110,11 @@ def get_coupled_variables(self, variables):
         variables.update(self._get_standard_overpotential_variables(eta_r))
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
 
+        if "main" in self.reaction:
+            variables.update(
+                self._get_standard_sei_film_overpotential_variables(eta_sei)
+            )
+
         if (
             "Negative electrode" + self.reaction_name + " interfacial current density"
             in variables
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 5b09b36d46..c93308a348 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -110,6 +110,8 @@ def __init__(self, model, param, disc, solution, operating_condition):
         ]
         self.eta_r_av = solution["X-averaged reaction overpotential [V]"]
 
+        self.eta_sei_av = solution["X-averaged sei film overpotential [V]"]
+
         self.eta_e_av = solution["X-averaged electrolyte overpotential [V]"]
         self.delta_phi_s_av = solution["X-averaged solid phase ohmic losses [V]"]
 
@@ -229,7 +231,8 @@ def test_consistent(self):
             self.ocv_av(self.t)
             + self.eta_r_av(self.t)
             + self.eta_e_av(self.t)
-            + self.delta_phi_s_av(self.t),
+            + self.delta_phi_s_av(self.t)
+            + self.eta_sei_av(self.t),
             decimal=2,
         )
 
@@ -530,13 +533,17 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.j_p_av = solution[
             "X-averaged positive electrode interfacial current density"
         ]
-        self.j_n_sei = solution["Negative electrode sei interfacial current density"]
-        self.j_p_sei = solution["Positive electrode sei interfacial current density"]
+        self.j_n_sei = solution[
+            "Negative electrode scaled sei interfacial current density"
+        ]
+        self.j_p_sei = solution[
+            "Positive electrode scaled sei interfacial current density"
+        ]
         self.j_n_sei_av = solution[
-            "X-averaged negative electrode sei interfacial current density"
+            "X-averaged negative electrode scaled sei interfacial current density"
         ]
         self.j_p_sei_av = solution[
-            "X-averaged positive electrode sei interfacial current density"
+            "X-averaged positive electrode scaled sei interfacial current density"
         ]
 
         self.j0_n = solution["Negative electrode exchange current density"]
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index cc59dbedba..33d1c8654d 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -7,113 +7,113 @@
 import unittest
 
 
-class TestSPM(unittest.TestCase):
-    def test_basic_processing(self):
-        options = {"thermal": "isothermal"}
-        model = pybamm.lithium_ion.SPM(options)
-        modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
-
-    def test_basic_processing_1plus1D(self):
-        options = {"current collector": "potential pair", "dimensionality": 1}
-
-        model = pybamm.lithium_ion.SPM(options)
-        var = pybamm.standard_spatial_vars
-        var_pts = {
-            var.x_n: 5,
-            var.x_s: 5,
-            var.x_p: 5,
-            var.r_n: 5,
-            var.r_p: 5,
-            var.y: 5,
-            var.z: 5,
-        }
-        modeltest = tests.StandardModelTest(model, var_pts=var_pts)
-        modeltest.test_all(skip_output_tests=True)
-
-    def test_basic_processing_2plus1D(self):
-        options = {"current collector": "potential pair", "dimensionality": 2}
-        model = pybamm.lithium_ion.SPM(options)
-        var = pybamm.standard_spatial_vars
-        var_pts = {
-            var.x_n: 5,
-            var.x_s: 5,
-            var.x_p: 5,
-            var.r_n: 5,
-            var.r_p: 5,
-            var.y: 5,
-            var.z: 5,
-        }
-        modeltest = tests.StandardModelTest(model, var_pts=var_pts)
-        modeltest.test_all(skip_output_tests=True)
-
-    def test_optimisations(self):
-        options = {"thermal": "isothermal"}
-        model = pybamm.lithium_ion.SPM(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)
-        simp_and_python = optimtest.evaluate_model(simplify=True, to_python=True)
-        np.testing.assert_array_almost_equal(original, simplified)
-        np.testing.assert_array_almost_equal(original, using_known_evals)
-        np.testing.assert_array_almost_equal(original, simp_and_known)
-        np.testing.assert_array_almost_equal(original, simp_and_python)
-
-    def test_set_up(self):
-        model = pybamm.lithium_ion.SPM()
-        optimtest = tests.OptimisationsTest(model)
-        optimtest.set_up_model(simplify=False, to_python=True)
-        optimtest.set_up_model(simplify=True, to_python=True)
-        optimtest.set_up_model(simplify=False, to_python=False)
-        optimtest.set_up_model(simplify=True, to_python=False)
-
-    def test_charge(self):
-        options = {"thermal": "isothermal"}
-        model = pybamm.lithium_ion.SPM(options)
-        parameter_values = model.default_parameter_values
-        parameter_values.update({"Current function [A]": -1})
-        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
-        modeltest.test_all()
-
-    def test_zero_current(self):
-        options = {"thermal": "isothermal"}
-        model = pybamm.lithium_ion.SPM(options)
-        parameter_values = model.default_parameter_values
-        parameter_values.update({"Current function [A]": 0})
-        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
-        modeltest.test_all()
-
-    def test_thermal(self):
-        options = {"thermal": "lumped"}
-        model = pybamm.lithium_ion.SPM(options)
-        modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
-
-        options = {"thermal": "x-full"}
-        model = pybamm.lithium_ion.SPM(options)
-        modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
-
-    def test_particle_fast_diffusion(self):
-        options = {"particle": "fast diffusion"}
-        model = pybamm.lithium_ion.SPM(options)
-        modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
-
-    def test_surface_form_differential(self):
-        options = {"surface form": "differential"}
-        model = pybamm.lithium_ion.SPM(options)
-        modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
-
-    def test_surface_form_algebraic(self):
-        options = {"surface form": "algebraic"}
-        model = pybamm.lithium_ion.SPM(options)
-        modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
+# class TestSPM(unittest.TestCase):
+#     def test_basic_processing(self):
+#         options = {"thermal": "isothermal"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         modeltest = tests.StandardModelTest(model)
+#         modeltest.test_all()
+
+#     def test_basic_processing_1plus1D(self):
+#         options = {"current collector": "potential pair", "dimensionality": 1}
+
+#         model = pybamm.lithium_ion.SPM(options)
+#         var = pybamm.standard_spatial_vars
+#         var_pts = {
+#             var.x_n: 5,
+#             var.x_s: 5,
+#             var.x_p: 5,
+#             var.r_n: 5,
+#             var.r_p: 5,
+#             var.y: 5,
+#             var.z: 5,
+#         }
+#         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
+#         modeltest.test_all(skip_output_tests=True)
+
+#     def test_basic_processing_2plus1D(self):
+#         options = {"current collector": "potential pair", "dimensionality": 2}
+#         model = pybamm.lithium_ion.SPM(options)
+#         var = pybamm.standard_spatial_vars
+#         var_pts = {
+#             var.x_n: 5,
+#             var.x_s: 5,
+#             var.x_p: 5,
+#             var.r_n: 5,
+#             var.r_p: 5,
+#             var.y: 5,
+#             var.z: 5,
+#         }
+#         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
+#         modeltest.test_all(skip_output_tests=True)
+
+#     def test_optimisations(self):
+#         options = {"thermal": "isothermal"}
+#         model = pybamm.lithium_ion.SPM(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)
+#         simp_and_python = optimtest.evaluate_model(simplify=True, to_python=True)
+#         np.testing.assert_array_almost_equal(original, simplified)
+#         np.testing.assert_array_almost_equal(original, using_known_evals)
+#         np.testing.assert_array_almost_equal(original, simp_and_known)
+#         np.testing.assert_array_almost_equal(original, simp_and_python)
+
+#     def test_set_up(self):
+#         model = pybamm.lithium_ion.SPM()
+#         optimtest = tests.OptimisationsTest(model)
+#         optimtest.set_up_model(simplify=False, to_python=True)
+#         optimtest.set_up_model(simplify=True, to_python=True)
+#         optimtest.set_up_model(simplify=False, to_python=False)
+#         optimtest.set_up_model(simplify=True, to_python=False)
+
+#     def test_charge(self):
+#         options = {"thermal": "isothermal"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         parameter_values = model.default_parameter_values
+#         parameter_values.update({"Current function [A]": -1})
+#         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+#         modeltest.test_all()
+
+#     def test_zero_current(self):
+#         options = {"thermal": "isothermal"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         parameter_values = model.default_parameter_values
+#         parameter_values.update({"Current function [A]": 0})
+#         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+#         modeltest.test_all()
+
+#     def test_thermal(self):
+#         options = {"thermal": "lumped"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         modeltest = tests.StandardModelTest(model)
+#         modeltest.test_all()
+
+#         options = {"thermal": "x-full"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         modeltest = tests.StandardModelTest(model)
+#         modeltest.test_all()
+
+#     def test_particle_fast_diffusion(self):
+#         options = {"particle": "fast diffusion"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         modeltest = tests.StandardModelTest(model)
+#         modeltest.test_all()
+
+#     def test_surface_form_differential(self):
+#         options = {"surface form": "differential"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         modeltest = tests.StandardModelTest(model)
+#         modeltest.test_all()
+
+#     def test_surface_form_algebraic(self):
+#         options = {"surface form": "algebraic"}
+#         model = pybamm.lithium_ion.SPM(options)
+#         modeltest = tests.StandardModelTest(model)
+#         modeltest.test_all()
 
 
 class TestSPMWithSEI(unittest.TestCase):

From 7c2601bf39078be4293353f8a32b6ed6cfa804b8 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 14 May 2020 15:31:37 -0400
Subject: [PATCH 11/15] #984 flake8 and test

---
 .../test_lithium_ion/test_spm.py              | 214 +++++++++---------
 .../test_leading_order_conductivity.py        |   3 +-
 2 files changed, 108 insertions(+), 109 deletions(-)

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 33d1c8654d..cc59dbedba 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -7,113 +7,113 @@
 import unittest
 
 
-# class TestSPM(unittest.TestCase):
-#     def test_basic_processing(self):
-#         options = {"thermal": "isothermal"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         modeltest = tests.StandardModelTest(model)
-#         modeltest.test_all()
-
-#     def test_basic_processing_1plus1D(self):
-#         options = {"current collector": "potential pair", "dimensionality": 1}
-
-#         model = pybamm.lithium_ion.SPM(options)
-#         var = pybamm.standard_spatial_vars
-#         var_pts = {
-#             var.x_n: 5,
-#             var.x_s: 5,
-#             var.x_p: 5,
-#             var.r_n: 5,
-#             var.r_p: 5,
-#             var.y: 5,
-#             var.z: 5,
-#         }
-#         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
-#         modeltest.test_all(skip_output_tests=True)
-
-#     def test_basic_processing_2plus1D(self):
-#         options = {"current collector": "potential pair", "dimensionality": 2}
-#         model = pybamm.lithium_ion.SPM(options)
-#         var = pybamm.standard_spatial_vars
-#         var_pts = {
-#             var.x_n: 5,
-#             var.x_s: 5,
-#             var.x_p: 5,
-#             var.r_n: 5,
-#             var.r_p: 5,
-#             var.y: 5,
-#             var.z: 5,
-#         }
-#         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
-#         modeltest.test_all(skip_output_tests=True)
-
-#     def test_optimisations(self):
-#         options = {"thermal": "isothermal"}
-#         model = pybamm.lithium_ion.SPM(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)
-#         simp_and_python = optimtest.evaluate_model(simplify=True, to_python=True)
-#         np.testing.assert_array_almost_equal(original, simplified)
-#         np.testing.assert_array_almost_equal(original, using_known_evals)
-#         np.testing.assert_array_almost_equal(original, simp_and_known)
-#         np.testing.assert_array_almost_equal(original, simp_and_python)
-
-#     def test_set_up(self):
-#         model = pybamm.lithium_ion.SPM()
-#         optimtest = tests.OptimisationsTest(model)
-#         optimtest.set_up_model(simplify=False, to_python=True)
-#         optimtest.set_up_model(simplify=True, to_python=True)
-#         optimtest.set_up_model(simplify=False, to_python=False)
-#         optimtest.set_up_model(simplify=True, to_python=False)
-
-#     def test_charge(self):
-#         options = {"thermal": "isothermal"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         parameter_values = model.default_parameter_values
-#         parameter_values.update({"Current function [A]": -1})
-#         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
-#         modeltest.test_all()
-
-#     def test_zero_current(self):
-#         options = {"thermal": "isothermal"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         parameter_values = model.default_parameter_values
-#         parameter_values.update({"Current function [A]": 0})
-#         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
-#         modeltest.test_all()
-
-#     def test_thermal(self):
-#         options = {"thermal": "lumped"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         modeltest = tests.StandardModelTest(model)
-#         modeltest.test_all()
-
-#         options = {"thermal": "x-full"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         modeltest = tests.StandardModelTest(model)
-#         modeltest.test_all()
-
-#     def test_particle_fast_diffusion(self):
-#         options = {"particle": "fast diffusion"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         modeltest = tests.StandardModelTest(model)
-#         modeltest.test_all()
-
-#     def test_surface_form_differential(self):
-#         options = {"surface form": "differential"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         modeltest = tests.StandardModelTest(model)
-#         modeltest.test_all()
-
-#     def test_surface_form_algebraic(self):
-#         options = {"surface form": "algebraic"}
-#         model = pybamm.lithium_ion.SPM(options)
-#         modeltest = tests.StandardModelTest(model)
-#         modeltest.test_all()
+class TestSPM(unittest.TestCase):
+    def test_basic_processing(self):
+        options = {"thermal": "isothermal"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_basic_processing_1plus1D(self):
+        options = {"current collector": "potential pair", "dimensionality": 1}
+
+        model = pybamm.lithium_ion.SPM(options)
+        var = pybamm.standard_spatial_vars
+        var_pts = {
+            var.x_n: 5,
+            var.x_s: 5,
+            var.x_p: 5,
+            var.r_n: 5,
+            var.r_p: 5,
+            var.y: 5,
+            var.z: 5,
+        }
+        modeltest = tests.StandardModelTest(model, var_pts=var_pts)
+        modeltest.test_all(skip_output_tests=True)
+
+    def test_basic_processing_2plus1D(self):
+        options = {"current collector": "potential pair", "dimensionality": 2}
+        model = pybamm.lithium_ion.SPM(options)
+        var = pybamm.standard_spatial_vars
+        var_pts = {
+            var.x_n: 5,
+            var.x_s: 5,
+            var.x_p: 5,
+            var.r_n: 5,
+            var.r_p: 5,
+            var.y: 5,
+            var.z: 5,
+        }
+        modeltest = tests.StandardModelTest(model, var_pts=var_pts)
+        modeltest.test_all(skip_output_tests=True)
+
+    def test_optimisations(self):
+        options = {"thermal": "isothermal"}
+        model = pybamm.lithium_ion.SPM(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)
+        simp_and_python = optimtest.evaluate_model(simplify=True, to_python=True)
+        np.testing.assert_array_almost_equal(original, simplified)
+        np.testing.assert_array_almost_equal(original, using_known_evals)
+        np.testing.assert_array_almost_equal(original, simp_and_known)
+        np.testing.assert_array_almost_equal(original, simp_and_python)
+
+    def test_set_up(self):
+        model = pybamm.lithium_ion.SPM()
+        optimtest = tests.OptimisationsTest(model)
+        optimtest.set_up_model(simplify=False, to_python=True)
+        optimtest.set_up_model(simplify=True, to_python=True)
+        optimtest.set_up_model(simplify=False, to_python=False)
+        optimtest.set_up_model(simplify=True, to_python=False)
+
+    def test_charge(self):
+        options = {"thermal": "isothermal"}
+        model = pybamm.lithium_ion.SPM(options)
+        parameter_values = model.default_parameter_values
+        parameter_values.update({"Current function [A]": -1})
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all()
+
+    def test_zero_current(self):
+        options = {"thermal": "isothermal"}
+        model = pybamm.lithium_ion.SPM(options)
+        parameter_values = model.default_parameter_values
+        parameter_values.update({"Current function [A]": 0})
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all()
+
+    def test_thermal(self):
+        options = {"thermal": "lumped"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+        options = {"thermal": "x-full"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_particle_fast_diffusion(self):
+        options = {"particle": "fast diffusion"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
 
 
 class TestSPMWithSEI(unittest.TestCase):
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_leading_order_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_leading_order_conductivity.py
index 98650b6dac..7d27313b2c 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_leading_order_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_leading_order_conductivity.py
@@ -14,8 +14,7 @@ def test_public_functions(self):
         variables = {
             "Current collector current density": a,
             "X-averaged negative electrode potential": a,
-            "X-averaged negative electrode open circuit potential": a,
-            "X-averaged negative electrode reaction overpotential": a,
+            "X-averaged negative electrode surface potential difference": a,
         }
         submodel = pybamm.electrolyte_conductivity.LeadingOrder(param)
         std_tests = tests.StandardSubModelTests(submodel, variables)

From 3afb1e4663d682a7b4f674c35af02c1bae5b67e9 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Fri, 15 May 2020 11:50:39 -0400
Subject: [PATCH 12/15] #984 use total current in sei film overpotential

---
 examples/scripts/calendar_ageing.py           |  5 --
 .../interface/kinetics/base_kinetics.py       | 57 ++++++++++---------
 2 files changed, 29 insertions(+), 33 deletions(-)

diff --git a/examples/scripts/calendar_ageing.py b/examples/scripts/calendar_ageing.py
index b0bdf731b7..6d8aa529e8 100644
--- a/examples/scripts/calendar_ageing.py
+++ b/examples/scripts/calendar_ageing.py
@@ -51,12 +51,7 @@
             "X-averaged negative electrode sei interfacial current density [A.m-2]",
             "X-averaged negative electrode interfacial current density [A.m-2]",
         ],
-        [
-            "X-averaged negative electrode sei interfacial current density",
-            "X-averaged negative electrode interfacial current density",
-        ],
         "Sum of x-averaged negative electrode interfacial current densities",
-        "Sum of negative electrode interfacial current densities",
         "X-averaged electrolyte concentration",
     ],
 )
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index 19b41ba60c..09de1a958c 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -37,16 +37,17 @@ def get_fundamental_variables(self):
             and "main" in self.reaction
         ):
             j = pybamm.Variable(
-                self.domain + " electrode interfacial current density",
+                "Total "
+                + self.domain.lower()
+                + " electrode interfacial current density",
                 domain=self.domain.lower() + " electrode",
                 auxiliary_domains={"secondary": "current collector"},
             )
 
             variables = {
-                self.domain
-                + " electrode"
-                + self.reaction_name
-                + " interfacial current density variable": j
+                "Total "
+                + self.domain.lower()
+                + " electrode interfacial current density variable": j
             }
             return variables
         else:
@@ -76,10 +77,12 @@ def get_coupled_variables(self, variables):
             L_sei = variables[
                 "Total " + self.domain.lower() + " electrode sei thickness"
             ]
-            j = variables[
-                self.domain + " electrode interfacial current density variable"
+            j_tot = variables[
+                "Total "
+                + self.domain.lower()
+                + " electrode interfacial current density variable"
             ]
-            eta_sei = -j * L_sei * pybamm.sei_parameters.R_sei
+            eta_sei = -j_tot * L_sei * pybamm.sei_parameters.R_sei
         elif self.options["sei film resistance"] == "average":
             L_sei = variables[
                 "Total " + self.domain.lower() + " electrode sei thickness"
@@ -138,20 +141,19 @@ def set_algebraic(self, variables):
             self.options["sei film resistance"] == "distributed"
             and "main" in self.reaction
         ):
-            j_var = variables[
-                self.domain
-                + " electrode"
-                + self.reaction_name
-                + " interfacial current density variable"
+            j_tot_var = variables[
+                "Total "
+                + self.domain.lower()
+                + " electrode interfacial current density variable"
             ]
-            j = variables[
-                self.domain
-                + " electrode"
-                + self.reaction_name
-                + " interfacial current density"
+            j_tot = variables[
+                "Sum of "
+                + self.domain.lower()
+                + " electrode interfacial current densities"
             ]
-            # Algebraic equation to set the variable j_var equal to the reaction term j
-            self.algebraic[j_var] = j_var - j
+            # Algebraic equation to set the variable j_tot_var
+            # equal to the sum of currents j_tot
+            self.algebraic[j_tot_var] = j_tot_var - j_tot
 
     def set_initial_conditions(self, variables):
         if (
@@ -159,18 +161,17 @@ def set_initial_conditions(self, variables):
             and "main" in self.reaction
         ):
             param = self.param
-            j_var = variables[
-                self.domain
-                + " electrode"
-                + self.reaction_name
-                + " interfacial current density variable"
+            j_tot_var = variables[
+                "Total "
+                + self.domain.lower()
+                + " electrode interfacial current density variable"
             ]
             if self.domain == "Negative":
-                j_av_init = param.current_with_time / param.l_n
+                j_tot_av_init = param.current_with_time / param.l_n
             elif self.domain == "Positive":
-                j_av_init = -param.current_with_time / param.l_p
+                j_tot_av_init = -param.current_with_time / param.l_p
 
-            self.initial_conditions[j_var] = j_av_init
+            self.initial_conditions[j_tot_var] = j_tot_av_init
 
     def _get_dj_dc(self, variables):
         """

From 01d04529ca9a8fac39df16a604bd397939efebe1 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Fri, 15 May 2020 14:48:25 -0400
Subject: [PATCH 13/15] #984 fix parameters

---
 pybamm/parameters/parameter_sets.py           |  7 +++++--
 pybamm/parameters/parameter_values.py         |  9 ++++++---
 .../test_parameters/test_parameter_values.py  | 19 ++++++++++++++++++-
 3 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/pybamm/parameters/parameter_sets.py b/pybamm/parameters/parameter_sets.py
index 6fb6536c3b..b487f79c97 100644
--- a/pybamm/parameters/parameter_sets.py
+++ b/pybamm/parameters/parameter_sets.py
@@ -18,7 +18,7 @@
     "cathode": "nca_Kim2011",
     "electrolyte": "lipf6_Kim2011",
     "experiment": "1C_discharge_from_full_Kim2011",
-    "sei": "Example",
+    "sei": "example",
     "citation": "kim2011multi",
 }
 
@@ -30,6 +30,7 @@
     "cathode": "LiNiCoO2_Ecker2015",
     "electrolyte": "lipf6_Ecker2015",
     "experiment": "1C_discharge_from_full_Ecker2015",
+    "sei": "example",
     "citation": ["ecker2015i", "ecker2015ii", "richardson2020"],
 }
 
@@ -41,7 +42,7 @@
     "cathode": "lico2_Marquis2019",
     "electrolyte": "lipf6_Marquis2019",
     "experiment": "1C_discharge_from_full_Marquis2019",
-    "sei": "Example",
+    "sei": "example",
     "citation": "marquis2019asymptotic",
 }
 
@@ -53,6 +54,7 @@
     "cathode": "nmc_Chen2020",
     "electrolyte": "lipf6_Nyman2008",
     "experiment": "1C_discharge_from_full_Chen2020",
+    "sei": "example",
     "citation": "Chen2020",
 }
 
@@ -64,6 +66,7 @@
     "cathode": "NMC_UMBL_Mohtat2020",
     "electrolyte": "LiPF6_Mohtat2020",
     "experiment": "1C_charge_from_empty_Mohtat2020",
+    "sei": "example",
     "citation": "Mohtat2020",
 }
 #
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 5a6b992fd9..5c782ec341 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -487,9 +487,12 @@ def _process_symbol(self, symbol):
                 # Scalar inherits name (for updating parameters) and domain (for
                 # Broadcast)
                 return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
-            elif isinstance(value, pybamm.InputParameter):
-                value.domain = symbol.domain
-                return value
+            elif isinstance(value, pybamm.Symbol):
+                new_value = self.process_symbol(value)
+                new_value.domain = symbol.domain
+                return new_value
+            else:
+                raise TypeError("Cannot process parameter '{}'".format(value))
 
         elif isinstance(symbol, pybamm.FunctionParameter):
             new_children = [self.process_symbol(child) for child in symbol.children]
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index aa35fdecbd..e37855279b 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -235,8 +235,20 @@ def test_process_symbol(self):
             x = pybamm.Parameter("x")
             parameter_values.process_symbol(x)
 
+    def test_process_parameter_in_parameter(self):
+        parameter_values = pybamm.ParameterValues(
+            {"a": 2, "2a": pybamm.Parameter("a") * 2}
+        )
+
+        # process 2a parameter
+        a = pybamm.Parameter("2a")
+        processed_a = parameter_values.process_symbol(a)
+        self.assertEqual(processed_a.evaluate(), 4)
+
     def test_process_input_parameter(self):
-        parameter_values = pybamm.ParameterValues({"a": "[input]", "b": 3})
+        parameter_values = pybamm.ParameterValues(
+            {"a": "[input]", "b": 3, "c times 2": pybamm.InputParameter("c") * 2}
+        )
         # process input parameter
         a = pybamm.Parameter("a")
         processed_a = parameter_values.process_symbol(a)
@@ -252,6 +264,11 @@ def test_process_input_parameter(self):
         self.assertIsInstance(processed_add.children[1], pybamm.Scalar)
         self.assertEqual(processed_add.evaluate(inputs={"a": 4}), 7)
 
+        # process complex input parameter
+        c = pybamm.Parameter("c times 2")
+        processed_c = parameter_values.process_symbol(c)
+        self.assertEqual(processed_c.evaluate(inputs={"c": 5}), 10)
+
     def test_process_function_parameter(self):
         parameter_values = pybamm.ParameterValues(
             {

From 4f44cd3d95cbc37e1a3a4827c47da84ed5e84af6 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Fri, 15 May 2020 15:18:07 -0400
Subject: [PATCH 14/15] #984 fix parameter path??

---
 .../parameters/lithium-ion/seis/{Example => example}/README.md    | 0
 .../lithium-ion/seis/{Example => example}/parameters.csv          | 0
 2 files changed, 0 insertions(+), 0 deletions(-)
 rename pybamm/input/parameters/lithium-ion/seis/{Example => example}/README.md (100%)
 rename pybamm/input/parameters/lithium-ion/seis/{Example => example}/parameters.csv (100%)

diff --git a/pybamm/input/parameters/lithium-ion/seis/Example/README.md b/pybamm/input/parameters/lithium-ion/seis/example/README.md
similarity index 100%
rename from pybamm/input/parameters/lithium-ion/seis/Example/README.md
rename to pybamm/input/parameters/lithium-ion/seis/example/README.md
diff --git a/pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv b/pybamm/input/parameters/lithium-ion/seis/example/parameters.csv
similarity index 100%
rename from pybamm/input/parameters/lithium-ion/seis/Example/parameters.csv
rename to pybamm/input/parameters/lithium-ion/seis/example/parameters.csv

From 3f1629c68e8b8d918f286eb37081ab9849b8ef81 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Tue, 19 May 2020 12:41:53 -0400
Subject: [PATCH 15/15] #984 change sei nondimensionalization

---
 .../submodels/interface/base_interface.py     |  8 ++---
 .../inverse_kinetics/inverse_butler_volmer.py |  4 +--
 .../submodels/interface/sei/base_sei.py       | 28 +++++----------
 .../sei/electron_migration_limited.py         | 10 ++++--
 .../sei/interstitial_diffusion_limited.py     | 14 +++++---
 .../interface/sei/reaction_limited.py         | 10 ++++--
 .../sei/solvent_diffusion_limited.py          | 11 ++++--
 pybamm/parameters/sei_parameters.py           | 35 ++++++++++---------
 .../test_models/standard_output_tests.py      | 12 +++----
 9 files changed, 71 insertions(+), 61 deletions(-)

diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index cb84e63b98..ad8bc0420f 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -35,9 +35,9 @@ def __init__(self, param, domain, reaction):
         elif reaction == "lithium-ion oxygen":
             self.reaction_name = " oxygen"
             self.Reaction_icd = "Oxygen interfacial current density"
-        elif reaction == "scaled sei":
-            self.reaction_name = " scaled sei"
-            self.Reaction_icd = "Scaled sei interfacial current density"
+        elif reaction == "sei":
+            self.reaction_name = " sei"
+            self.Reaction_icd = "Sei interfacial current density"
         self.reaction = reaction
 
     def _get_exchange_current_density(self, variables):
@@ -175,7 +175,7 @@ def _get_number_of_electrons_in_reaction(self):
 
     def _get_electrolyte_reaction_signed_stoichiometry(self):
         "Returns the number of electrons in the reaction"
-        if self.reaction in ["lithium-ion main", "scaled sei"]:
+        if self.reaction in ["lithium-ion main", "sei"]:
             # Both the main reaction current contribute to the electrolyte reaction
             # current
             return pybamm.Scalar(1), pybamm.Scalar(1)
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 6531adfe24..110ab39bba 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -131,9 +131,7 @@ def get_coupled_variables(self, variables):
             + self.domain.lower()
             + " electrode total interfacial current density"
         ]
-        j_sei = variables[
-            self.domain + " electrode scaled sei interfacial current density"
-        ]
+        j_sei = variables[self.domain + " electrode sei interfacial current density"]
         j = j_tot - j_sei
 
         variables.update(self._get_standard_interfacial_current_variables(j))
diff --git a/pybamm/models/submodels/interface/sei/base_sei.py b/pybamm/models/submodels/interface/sei/base_sei.py
index 2feba30c61..d457344098 100644
--- a/pybamm/models/submodels/interface/sei/base_sei.py
+++ b/pybamm/models/submodels/interface/sei/base_sei.py
@@ -19,7 +19,11 @@ class BaseModel(BaseInterface):
     """
 
     def __init__(self, param, domain):
-        reaction = "scaled sei"
+        if domain == "Positive" and not isinstance(self, pybamm.sei.NoSEI):
+            raise NotImplementedError(
+                "SEI models are not implemented for the positive electrode"
+            )
+        reaction = "sei"
         super().__init__(param, domain, reaction)
 
     def _get_standard_thickness_variables(self, L_inner, L_outer):
@@ -117,20 +121,10 @@ def _get_standard_reaction_variables(self, j_inner, j_outer):
             variables : dict
                 The variables which can be derived from the SEI thicknesses.
         """
-        # Set scales to one for the "no SEI" model so that they are not required
-        # by parameter values in general
-        if isinstance(self, pybamm.sei.NoSEI):
-            j_scale = 1
-            Gamma_SEI = 1
-        else:
-            sp = pybamm.sei_parameters
-            j_scale = (
-                sp.F * sp.L_sei_0_dim / sp.V_bar_inner_dimensional / sp.tau_discharge
-            )
-            if self.domain == "Negative":
-                Gamma_SEI = sp.Gamma_SEI_n
-            elif self.domain == "Positive":
-                Gamma_SEI = sp.Gamma_SEI_p
+        if self.domain == "Negative":
+            j_scale = self.param.interfacial_current_scale_n
+        elif self.domain == "Positive":
+            j_scale = self.param.interfacial_current_scale_p
         j_i_av = pybamm.x_average(j_inner)
         j_o_av = pybamm.x_average(j_outer)
 
@@ -163,10 +157,6 @@ def _get_standard_reaction_variables(self, j_inner, j_outer):
             "X-averaged "
             + domain
             + " sei interfacial current density [A.m-2]": j_sei_av * j_scale,
-            Domain + " scaled sei interfacial current density": j_sei * Gamma_SEI,
-            "X-averaged "
-            + domain
-            + " scaled sei interfacial current density": j_sei_av * Gamma_SEI,
         }
 
         return variables
diff --git a/pybamm/models/submodels/interface/sei/electron_migration_limited.py b/pybamm/models/submodels/interface/sei/electron_migration_limited.py
index 9d76dd80e3..f351666311 100644
--- a/pybamm/models/submodels/interface/sei/electron_migration_limited.py
+++ b/pybamm/models/submodels/interface/sei/electron_migration_limited.py
@@ -35,8 +35,9 @@ def get_coupled_variables(self, variables):
         ]
         phi_s_n = variables[self.domain + " electrode potential"]
 
-        C_sei = pybamm.sei_parameters.C_sei_electron
         U_inner = pybamm.sei_parameters.U_inner_electron
+        if self.domain == "Negative":
+            C_sei = pybamm.sei_parameters.C_sei_electron_n
 
         j_sei = (phi_s_n - U_inner) / (C_sei * L_sei_inner)
 
@@ -66,8 +67,13 @@ def set_rhs(self, variables):
         j_outer = variables["Outer " + domain + " sei interfacial current density"]
 
         v_bar = pybamm.sei_parameters.v_bar
+        if self.domain == "Negative":
+            Gamma_SEI = pybamm.sei_parameters.Gamma_SEI_n
 
-        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
+        self.rhs = {
+            L_inner: -Gamma_SEI * j_inner,
+            L_outer: -v_bar * Gamma_SEI * j_outer,
+        }
 
     def set_initial_conditions(self, variables):
         domain = self.domain.lower() + " electrode"
diff --git a/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py b/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
index 9e15caa522..b962e465cb 100644
--- a/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/sei/interstitial_diffusion_limited.py
@@ -34,11 +34,11 @@ def get_coupled_variables(self, variables):
             "Inner " + self.domain.lower() + " electrode sei thickness"
         ]
         phi_s_n = variables[self.domain + " electrode potential"]
-        phi_e_n = variables[self.domain + " electrolyte potential"]
 
-        C_sei = pybamm.sei_parameters.C_sei_inter
+        if self.domain == "Negative":
+            C_sei = pybamm.sei_parameters.C_sei_inter_n
 
-        j_sei = -pybamm.exp(-(phi_s_n - phi_e_n)) / (C_sei * L_sei_inner)
+        j_sei = -pybamm.exp(-phi_s_n) / (C_sei * L_sei_inner)
 
         alpha = 0.5
         j_inner = alpha * j_sei
@@ -67,7 +67,13 @@ def set_rhs(self, variables):
 
         v_bar = pybamm.sei_parameters.v_bar
 
-        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
+        if self.domain == "Negative":
+            Gamma_SEI = pybamm.sei_parameters.Gamma_SEI_n
+
+        self.rhs = {
+            L_inner: -Gamma_SEI * j_inner,
+            L_outer: -v_bar * Gamma_SEI * j_outer,
+        }
 
     def set_initial_conditions(self, variables):
         domain = self.domain.lower() + " electrode"
diff --git a/pybamm/models/submodels/interface/sei/reaction_limited.py b/pybamm/models/submodels/interface/sei/reaction_limited.py
index 69aa589e8c..f312644cd5 100644
--- a/pybamm/models/submodels/interface/sei/reaction_limited.py
+++ b/pybamm/models/submodels/interface/sei/reaction_limited.py
@@ -48,10 +48,11 @@ def get_coupled_variables(self, variables):
             ]
         L_sei = variables["Total " + self.domain.lower() + " electrode sei thickness"]
 
-        C_sei = pybamm.sei_parameters.C_sei_reaction
         R_sei = pybamm.sei_parameters.R_sei
         alpha = 0.5
         # alpha = pybamm.sei_parameters.alpha
+        if self.domain == "Negative":
+            C_sei = pybamm.sei_parameters.C_sei_reaction_n
 
         # need to revise for thermal case
         j_sei = -(1 / C_sei) * pybamm.exp(
@@ -83,8 +84,13 @@ def set_rhs(self, variables):
         j_outer = variables["Outer " + domain + " sei interfacial current density"]
 
         v_bar = pybamm.sei_parameters.v_bar
+        if self.domain == "Negative":
+            Gamma_SEI = pybamm.sei_parameters.Gamma_SEI_n
 
-        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
+        self.rhs = {
+            L_inner: -Gamma_SEI * j_inner,
+            L_outer: -v_bar * Gamma_SEI * j_outer,
+        }
 
     def set_initial_conditions(self, variables):
         domain = self.domain.lower() + " electrode"
diff --git a/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py b/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
index a1e638ba34..fc2c3f8b8a 100644
--- a/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/sei/solvent_diffusion_limited.py
@@ -34,7 +34,8 @@ def get_coupled_variables(self, variables):
             "Outer " + self.domain.lower() + " electrode sei thickness"
         ]
 
-        C_sei = pybamm.sei_parameters.C_sei_solvent
+        if self.domain == "Negative":
+            C_sei = pybamm.sei_parameters.C_sei_solvent_n
 
         j_sei = -1 / (C_sei * L_sei_outer)
 
@@ -65,7 +66,13 @@ def set_rhs(self, variables):
 
         v_bar = pybamm.sei_parameters.v_bar
 
-        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
+        if self.domain == "Negative":
+            Gamma_SEI = pybamm.sei_parameters.Gamma_SEI_n
+
+        self.rhs = {
+            L_inner: -Gamma_SEI * j_inner,
+            L_outer: -v_bar * Gamma_SEI * j_outer,
+        }
 
     def set_initial_conditions(self, variables):
         domain = self.domain.lower() + " electrode"
diff --git a/pybamm/parameters/sei_parameters.py b/pybamm/parameters/sei_parameters.py
index ef1106f7b0..d4532f5757 100644
--- a/pybamm/parameters/sei_parameters.py
+++ b/pybamm/parameters/sei_parameters.py
@@ -51,24 +51,25 @@
 L_x = pybamm.standard_parameters_lithium_ion.L_x
 
 i_typ = pybamm.electrical_parameters.i_typ
+j_scale_n = pybamm.standard_parameters_lithium_ion.interfacial_current_scale_n
+j_scale_p = pybamm.standard_parameters_lithium_ion.interfacial_current_scale_p
 
-C_sei_reaction = (
-    F * L_sei_0_dim / (m_sei_dimensional * tau_discharge * V_bar_inner_dimensional)
-) * pybamm.exp(-(F * U_n_ref / (2 * R * T_ref)))
 
-C_sei_solvent = L_sei_0_dim ** 2 / (
-    c_sol_dimensional * V_bar_inner_dimensional * D_sol_dimensional * tau_discharge
+C_sei_reaction_n = (j_scale_n / m_sei_dimensional) * pybamm.exp(
+    -(F * U_n_ref / (2 * R * T_ref))
 )
-
-C_sei_electron = (
-    F ** 2
-    * L_sei_0_dim ** 2
-    / (kappa_inner_dimensional * V_bar_inner_dimensional * R * T_ref * tau_discharge)
+C_sei_reaction_p = (j_scale_p / m_sei_dimensional) * pybamm.exp(
+    -(F * U_n_ref / (2 * R * T_ref))
 )
 
-C_sei_inter = L_sei_0_dim ** 2 / (
-    D_li_dimensional * c_li_0_dimensional * V_bar_inner_dimensional * tau_discharge
-)
+C_sei_solvent_n = j_scale_n * L_sei_0_dim / (c_sol_dimensional * F * D_sol_dimensional)
+C_sei_solvent_p = j_scale_p * L_sei_0_dim / (c_sol_dimensional * F * D_sol_dimensional)
+
+C_sei_electron_n = j_scale_n * F * L_sei_0_dim / (kappa_inner_dimensional * R * T_ref)
+C_sei_electron_p = j_scale_p * F * L_sei_0_dim / (kappa_inner_dimensional * R * T_ref)
+
+C_sei_inter_n = j_scale_n * L_sei_0_dim / (D_li_dimensional * c_li_0_dimensional * F)
+C_sei_inter_p = j_scale_p * L_sei_0_dim / (D_li_dimensional * c_li_0_dimensional * F)
 
 U_inner_electron = F * U_inner_dimensional / R / T_ref
 
@@ -80,10 +81,10 @@
 L_outer_0 = L_outer_0_dim / L_sei_0_dim
 
 # ratio of SEI reaction scale to intercalation reaction
-Gamma_SEI_n = (F * L_sei_0_dim * a_n * L_x) / (
-    V_bar_inner_dimensional * i_typ * tau_discharge
+Gamma_SEI_n = (V_bar_inner_dimensional * i_typ * tau_discharge) / (
+    F * L_sei_0_dim * a_n * L_x
 )
 
-Gamma_SEI_p = (F * L_sei_0_dim * a_p * L_x) / (
-    V_bar_inner_dimensional * i_typ * tau_discharge
+Gamma_SEI_p = (V_bar_inner_dimensional * i_typ * tau_discharge) / (
+    F * L_sei_0_dim * a_p * L_x
 )
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index c93308a348..3732f01623 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -533,17 +533,13 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.j_p_av = solution[
             "X-averaged positive electrode interfacial current density"
         ]
-        self.j_n_sei = solution[
-            "Negative electrode scaled sei interfacial current density"
-        ]
-        self.j_p_sei = solution[
-            "Positive electrode scaled sei interfacial current density"
-        ]
+        self.j_n_sei = solution["Negative electrode sei interfacial current density"]
+        self.j_p_sei = solution["Positive electrode sei interfacial current density"]
         self.j_n_sei_av = solution[
-            "X-averaged negative electrode scaled sei interfacial current density"
+            "X-averaged negative electrode sei interfacial current density"
         ]
         self.j_p_sei_av = solution[
-            "X-averaged positive electrode scaled sei interfacial current density"
+            "X-averaged positive electrode sei interfacial current density"
         ]
 
         self.j0_n = solution["Negative electrode exchange current density"]