diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 2a98f67f9e..99d0b5fa60 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -124,6 +124,11 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "SEI porosity change" : str
                 Whether to include porosity change due to SEI formation, can be "false"
                 (default) or "true".
+            * "stress induced diffusion" : str
+                Whether to includes stress induced diffusion, can be "false" or "true".
+                The default is "false" if "particle mechanics" is "none" and "true"
+                otherwise. A 2-tuple can be provided for different behaviour in negative
+                and positive electrodes.
             * "surface form" : str
                 Whether to use the surface formulation of the problem. Can be "false"
                 (default), "differential" or "algebraic".
@@ -186,6 +191,7 @@ def __init__(self, extra_options):
             ],
             "SEI film resistance": ["none", "distributed", "average"],
             "SEI porosity change": ["true", "false"],
+            "stress induced diffusion": ["true", "false"],
             "surface form": ["false", "differential", "algebraic"],
             "thermal": ["isothermal", "lumped", "x-lumped", "x-full"],
             "total interfacial current density as a state": ["true", "false"],
@@ -238,7 +244,7 @@ def __init__(self, extra_options):
         # The "SEI film resistance" option will still be overridden by extra_options if
         # provided
 
-        # Change the default for swelling based on which LAM option is
+        # Change the default for particle mechanics based on which LAM option is
         # provided
         # return "none" if option not given
         lam_option = extra_options.get("loss of active material", "none")
@@ -246,9 +252,32 @@ def __init__(self, extra_options):
             default_options["particle mechanics"] = "swelling only"
         else:
             default_options["particle mechanics"] = "none"
-        # The "SEI film resistance" option will still be overridden by extra_options if
+        # The "particle mechanics" option will still be overridden by extra_options if
         # provided
 
+        # Change the default for stress induced diffusion based on which particle
+        # mechanics option is provided
+        mechanics_options = extra_options.get("particle mechanics", "none")
+        if mechanics_options == "none":
+            default_options["stress induced diffusion"] = "false"
+        else:
+            if isinstance(self.options["particle mechanics"], str):
+                mech_left = self.options["particle mechanics"]
+                mech_right = self.options["particle mechanics"]
+            else:
+                mech_left, mech_right = self.options["particle mechanics"]
+            if mech_left == "none":
+                default_left = "false"
+            else:
+                default_left = "true"
+            if mech_right == "none":
+                default_right = "false"
+            else:
+                default_right = "true"
+            default_options["stress induced diffusion"] = (default_left, default_right)
+        # The "stress induced diffusion" option will still be overridden by
+        # extra_options if provided
+
         options = pybamm.FuzzyDict(default_options)
         # any extra options overwrite the default options
         for name, opt in extra_options.items():
@@ -283,10 +312,6 @@ def __init__(self, extra_options):
 
         # Options not yet compatible with particle-size distributions
         if options["particle size"] == "distribution":
-            if options["SEI"] != "none":
-                raise NotImplementedError(
-                    "SEI submodels do not yet support particle-size distributions."
-                )
             if options["lithium plating"] != "none":
                 raise NotImplementedError(
                     "Lithium plating submodels do not yet support particle-size "
@@ -299,9 +324,18 @@ def __init__(self, extra_options):
                 )
             if options["particle shape"] != "spherical":
                 raise NotImplementedError(
-                    "Particle shape must be 'spherical' for particle-size distributions"
+                    "Particle shape must be 'spherical' for particle-size distribution"
                     " submodels."
                 )
+            if options["SEI"] != "none":
+                raise NotImplementedError(
+                    "SEI submodels do not yet support particle-size distributions."
+                )
+            if options["stress induced diffusion"] == "true":
+                raise NotImplementedError(
+                    "Stress induced diffusion cannot yet be included in "
+                    "particle-size distributions."
+                )
             if options["thermal"] == "x-full":
                 raise NotImplementedError(
                     "X-full thermal submodels do not yet support particle-size"
diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py
index 7d9b398270..514a41b474 100644
--- a/pybamm/models/standard_variables.py
+++ b/pybamm/models/standard_variables.py
@@ -175,12 +175,12 @@ def __init__(self):
             auxiliary_domains={"secondary": "current collector"},
             bounds=(0, 1),
         )
-        self.c_s_n_rxav = pybamm.Variable(
+        self.c_s_n_av = pybamm.Variable(
             "Average negative particle concentration",
             domain="current collector",
             bounds=(0, 1),
         )
-        self.c_s_p_rxav = pybamm.Variable(
+        self.c_s_p_av = pybamm.Variable(
             "Average positive particle concentration",
             domain="current collector",
             bounds=(0, 1),
@@ -220,11 +220,11 @@ def __init__(self):
             domain="positive electrode",
             auxiliary_domains={"secondary": "current collector"},
         )
-        self.q_s_n_rxav = pybamm.Variable(
+        self.q_s_n_av = pybamm.Variable(
             "Average negative particle concentration gradient",
             domain="current collector",
         )
-        self.q_s_p_rxav = pybamm.Variable(
+        self.q_s_p_av = pybamm.Variable(
             "Average positive particle concentration gradient",
             domain="current collector",
         )
diff --git a/pybamm/models/submodels/particle/base_fickian.py b/pybamm/models/submodels/particle/base_fickian.py
new file mode 100644
index 0000000000..c5341b5cb8
--- /dev/null
+++ b/pybamm/models/submodels/particle/base_fickian.py
@@ -0,0 +1,63 @@
+#
+# Base class for particles with Fickian diffusion
+#
+import pybamm
+from base_particle import BaseParticle
+
+
+class BaseFickian(BaseParticle):
+    """
+    Base class for molar conservation in particles with Fickian diffusion.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    domain : str
+        The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
+
+    **Extends:** :class:`pybamm.particle.BaseParticle`
+    """
+
+    def __init__(self, param, domain, options):
+        super().__init__(param, domain, options)
+
+        if self.options["stress induced diffusion"] == "true":
+            pybamm.citations.register("Ai2019")
+            pybamm.citations.register("Deshpande2012")
+
+    def _get_effective_diffusivity(self, c, T):
+        if self.domain == "Negative":
+            D = self.param.D_n(c, T)
+        elif self.domain == "Positive":
+            D = self.param.D_p(c, T)
+
+        # Account for stress induced diffusion
+        if self.options["stress induced diffusion"] == "true":
+            if self.domain == "Negative":
+                theta = self.param.theta_n
+                c_0 = self.param.c_0_n
+            elif self.domain == "Positive":
+                theta = self.param.theta_p
+                c_0 = self.param.c_0_p
+
+            D_eff = D * (1 + theta * (c - c_0) / (1 + self.param.Theta * T))
+        else:
+            D_eff = D
+
+        return D_eff
+
+    def _get_standard_flux_variables(self, N_s, N_s_xav, D_eff):
+        variables = {
+            self.domain + " particle flux": N_s,
+            "X-averaged " + self.domain.lower() + " particle flux": N_s_xav,
+            self.domain + " effective diffusivity": D_eff,
+            "X-averaged "
+            + self.domain.lower()
+            + " effective diffusivity": pybamm.x_average(D_eff),
+        }
+
+        return variables
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index cdbda45a31..f12df7af1a 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -14,12 +14,15 @@ class BaseParticle(pybamm.BaseSubModel):
         The parameters to use for this submodel
     domain : str
         The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, options=None):
+        super().__init__(param, domain, options=options)
 
     def _get_standard_concentration_variables(
         self, c_s, c_s_xav=None, c_s_rav=None, c_s_av=None, c_s_surf=None
@@ -148,4 +151,3 @@ def _get_standard_flux_variables(self, N_s, N_s_xav):
         }
 
         return variables
-
diff --git a/pybamm/models/submodels/particle/no_distribution/fickian_diffusion.py b/pybamm/models/submodels/particle/no_distribution/fickian_diffusion.py
index ff9ee7ca34..08a5dc7d43 100644
--- a/pybamm/models/submodels/particle/no_distribution/fickian_diffusion.py
+++ b/pybamm/models/submodels/particle/no_distribution/fickian_diffusion.py
@@ -2,10 +2,10 @@
 # Class for particles with Fickian diffusion and x-dependence
 #
 import pybamm
-from ..base_particle import BaseParticle
+from ..base_fickian import BaseFickian
 
 
-class FickianDiffusion(BaseParticle):
+class FickianDiffusion(BaseFickian):
     """
     Class for molar conservation in particles, employing Fick's law, and allowing
     variation in the electrode domain. I.e., the concentration varies with r
@@ -17,12 +17,15 @@ class FickianDiffusion(BaseParticle):
         The parameters to use for this submodel
     domain : str
         The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.particle.BaseParticle`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, options):
+        super().__init__(param, domain, options)
 
     def get_fundamental_variables(self):
         if self.domain == "Negative":
@@ -42,13 +45,10 @@ def get_coupled_variables(self, variables):
             [self.domain.lower() + " particle"],
         )
 
-        if self.domain == "Negative":
-            N_s = -self.param.D_n(c_s, T) * pybamm.grad(c_s)
-
-        elif self.domain == "Positive":
-            N_s = -self.param.D_p(c_s, T) * pybamm.grad(c_s)
+        D_eff = self._get_effective_diffusivity(c_s, T)
+        N_s = -D_eff * pybamm.grad(c_s)
 
-        variables.update(self._get_standard_flux_variables(N_s, N_s))
+        variables.update(self._get_standard_flux_variables(N_s, N_s, D_eff))
         variables.update(self._get_total_concentration_variables(variables))
 
         return variables
@@ -67,21 +67,12 @@ def set_rhs(self, variables):
     def set_boundary_conditions(self, variables):
 
         c_s = variables[self.domain + " particle concentration"]
-        T = pybamm.PrimaryBroadcast(
-            variables[self.domain + " electrode temperature"],
-            c_s.domain,
-        )
+        D_eff = variables[self.domain + " effective diffusivity"]
         j = variables[self.domain + " electrode interfacial current density"]
         R = variables[self.domain + " particle radius"]
 
         if self.domain == "Negative":
-            rbc = (
-                -self.param.C_n
-                * j
-                * R
-                / self.param.a_R_n
-                / pybamm.surf(self.param.D_n(c_s, T))
-            )
+            rbc = -self.param.C_n * j * R / self.param.a_R_n / pybamm.surf(D_eff)
 
         elif self.domain == "Positive":
             rbc = (
@@ -90,7 +81,7 @@ def set_boundary_conditions(self, variables):
                 * R
                 / self.param.a_R_p
                 / self.param.gamma_p
-                / pybamm.surf(self.param.D_p(c_s, T))
+                / pybamm.surf(D_eff)
             )
 
         self.boundary_conditions = {
diff --git a/pybamm/models/submodels/particle/no_distribution/polynomial_profile.py b/pybamm/models/submodels/particle/no_distribution/polynomial_profile.py
index 9217ee1e7c..ee5fa95cee 100644
--- a/pybamm/models/submodels/particle/no_distribution/polynomial_profile.py
+++ b/pybamm/models/submodels/particle/no_distribution/polynomial_profile.py
@@ -3,14 +3,14 @@
 #
 import pybamm
 
-from ..base_particle import BaseParticle
+from ..base_fickian import BaseFickian
 
 
-class PolynomialProfile(BaseParticle):
+class PolynomialProfile(BaseFickian):
     """
-    Class for molar conservation in particles, assuming a polynomial
-    concentration profile in r, and allowing variation in the electrode domain.
-    Model equations from [1]_.
+    Class for molar conservation in particles employing Fick's
+    law, assuming a polynomial concentration profile in r, and allowing variation
+    in the electrode domain. Model equations from [1]_.
 
     Parameters
     ----------
@@ -21,6 +21,9 @@ class PolynomialProfile(BaseParticle):
     name : str
         The name of the polynomial approximation to be used. Can be "uniform
         profile", "quadratic profile" or "quartic profile".
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     References
     ----------
@@ -31,8 +34,8 @@ class PolynomialProfile(BaseParticle):
     **Extends:** :class:`pybamm.particle.BaseParticle`
     """
 
-    def __init__(self, param, domain, name):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, name, options):
+        super().__init__(param, domain, options)
         self.name = name
 
         pybamm.citations.register("Subramanian2005")
@@ -133,6 +136,7 @@ def get_coupled_variables(self, variables):
             variables[self.domain + " electrode temperature"],
             [self.domain.lower() + " particle"],
         )
+        D_eff = self._get_effective_diffusivity(c_s, T)
 
         # Set flux depending on polynomial order
         if self.name == "uniform profile":
@@ -152,10 +156,10 @@ def get_coupled_variables(self, variables):
             # The flux may be computed directly from the polynomial for c
             if self.domain == "Negative":
                 r = pybamm.standard_spatial_vars.r_n
-                N_s = -self.param.D_n(c_s, T) * 5 * (c_s_surf - c_s_rav) * r
+                N_s = -D_eff * 5 * (c_s_surf - c_s_rav) * r
             elif self.domain == "Positive":
                 r = pybamm.standard_spatial_vars.r_p
-                N_s = -self.param.D_p(c_s, T) * 5 * (c_s_surf - c_s_rav) * r
+                N_s = -D_eff * 5 * (c_s_surf - c_s_rav) * r
             N_s_xav = pybamm.x_average(N_s)
         elif self.name == "quartic profile":
             q_s_rav = variables[
@@ -164,19 +168,19 @@ def get_coupled_variables(self, variables):
             # The flux may be computed directly from the polynomial for c
             if self.domain == "Negative":
                 r = pybamm.standard_spatial_vars.r_n
-                N_s = -self.param.D_n(c_s, T) * (
+                N_s = -D_eff * (
                     (-70 * c_s_surf + 20 * q_s_rav + 70 * c_s_rav) * r
                     + (105 * c_s_surf - 28 * q_s_rav - 105 * c_s_rav) * r ** 3
                 )
             elif self.domain == "Positive":
                 r = pybamm.standard_spatial_vars.r_p
-                N_s = -self.param.D_p(c_s, T) * (
+                N_s = -D_eff * (
                     (-70 * c_s_surf + 20 * q_s_rav + 70 * c_s_rav) * r
                     + (105 * c_s_surf - 28 * q_s_rav - 105 * c_s_rav) * r ** 3
                 )
             N_s_xav = pybamm.x_average(N_s)
 
-        variables.update(self._get_standard_flux_variables(N_s, N_s_xav))
+        variables.update(self._get_standard_flux_variables(N_s, N_s_xav, D_eff))
         variables.update(self._get_total_concentration_variables(variables))
 
         return variables
@@ -202,12 +206,13 @@ def set_rhs(self, variables):
             c_s_rav = variables[
                 "R-averaged " + self.domain.lower() + " particle concentration"
             ]
-            T = variables[self.domain + " electrode temperature"]
+            D_eff = variables[self.domain + " effective diffusivity"]
+
             if self.domain == "Negative":
                 self.rhs.update(
                     {
                         q_s_rav: -30
-                        * self.param.D_n(c_s_rav, T)
+                        * pybamm.r_average(D_eff)
                         * q_s_rav
                         / self.param.C_n
                         - 45 * j / self.param.a_R_n / 2
@@ -217,7 +222,7 @@ def set_rhs(self, variables):
                 self.rhs.update(
                     {
                         q_s_rav: -30
-                        * self.param.D_p(c_s_rav, T)
+                        * pybamm.r_average(D_eff)
                         * q_s_rav
                         / self.param.C_p
                         - 45 * j / self.param.a_R_p / self.param.gamma_p / 2
@@ -229,9 +234,10 @@ def set_algebraic(self, variables):
         c_s_rav = variables[
             "R-averaged " + self.domain.lower() + " particle concentration"
         ]
+        D_eff = variables[self.domain + " effective diffusivity"]
         j = variables[self.domain + " electrode interfacial current density"]
-        T = variables[self.domain + " electrode temperature"]
         R = variables[self.domain + " particle radius"]
+
         if self.name == "uniform profile":
             # No algebraic equations since we only solve for the average concentration
             pass
@@ -239,13 +245,13 @@ def set_algebraic(self, variables):
             # We solve an algebraic equation for the surface concentration
             if self.domain == "Negative":
                 self.algebraic = {
-                    c_s_surf: self.param.D_n(c_s_surf, T) * (c_s_surf - c_s_rav)
+                    c_s_surf: pybamm.surf(D_eff) * (c_s_surf - c_s_rav)
                     + self.param.C_n * (j * R / self.param.a_R_n / 5)
                 }
 
             elif self.domain == "Positive":
                 self.algebraic = {
-                    c_s_surf: self.param.D_p(c_s_surf, T) * (c_s_surf - c_s_rav)
+                    c_s_surf: pybamm.surf(D_eff) * (c_s_surf - c_s_rav)
                     + self.param.C_p
                     * (j * R / self.param.a_R_p / self.param.gamma_p / 5)
                 }
@@ -257,14 +263,14 @@ def set_algebraic(self, variables):
             ]
             if self.domain == "Negative":
                 self.algebraic = {
-                    c_s_surf: self.param.D_n(c_s_surf, T)
+                    c_s_surf: pybamm.surf(D_eff)
                     * (35 * (c_s_surf - c_s_rav) - 8 * q_s_rav)
                     + self.param.C_n * (j * R / self.param.a_R_n)
                 }
 
             elif self.domain == "Positive":
                 self.algebraic = {
-                    c_s_surf: self.param.D_p(c_s_surf, T)
+                    c_s_surf: pybamm.surf(D_eff)
                     * (35 * (c_s_surf - c_s_rav) - 8 * q_s_rav)
                     + self.param.C_p * (j * R / self.param.a_R_p / self.param.gamma_p)
                 }
diff --git a/pybamm/models/submodels/particle/no_distribution/x_averaged_fickian_diffusion.py b/pybamm/models/submodels/particle/no_distribution/x_averaged_fickian_diffusion.py
index ea4370af6d..173e5ac179 100644
--- a/pybamm/models/submodels/particle/no_distribution/x_averaged_fickian_diffusion.py
+++ b/pybamm/models/submodels/particle/no_distribution/x_averaged_fickian_diffusion.py
@@ -3,10 +3,10 @@
 #
 import pybamm
 
-from ..base_particle import BaseParticle
+from ..base_fickian import BaseFickian
 
 
-class XAveragedFickianDiffusion(BaseParticle):
+class XAveragedFickianDiffusion(BaseFickian):
     """
     Class for molar conservation in a single x-averaged particle, employing Fick's
     law. I.e., the concentration varies with r (internal spherical coordinate)
@@ -18,12 +18,15 @@ class XAveragedFickianDiffusion(BaseParticle):
         The parameters to use for this submodel
     domain : str
         The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.particle.BaseParticle`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, options):
+        super().__init__(param, domain, options)
 
     def get_fundamental_variables(self):
         if self.domain == "Negative":
@@ -47,15 +50,15 @@ def get_coupled_variables(self, variables):
             [self.domain.lower() + " particle"],
         )
 
-        if self.domain == "Negative":
-            N_s_xav = -self.param.D_n(c_s_xav, T_xav) * pybamm.grad(c_s_xav)
-
-        elif self.domain == "Positive":
-            N_s_xav = -self.param.D_p(c_s_xav, T_xav) * pybamm.grad(c_s_xav)
+        D_eff_xav = self._get_effective_diffusivity(c_s_xav, T_xav)
+        N_s_xav = -D_eff_xav * pybamm.grad(c_s_xav)
 
+        D_eff = pybamm.SecondaryBroadcast(
+            D_eff_xav, [self._domain.lower() + " electrode"]
+        )
         N_s = pybamm.SecondaryBroadcast(N_s_xav, [self._domain.lower() + " electrode"])
 
-        variables.update(self._get_standard_flux_variables(N_s, N_s_xav))
+        variables.update(self._get_standard_flux_variables(N_s, N_s_xav, D_eff))
         variables.update(self._get_total_concentration_variables(variables))
 
         return variables
@@ -76,10 +79,9 @@ def set_boundary_conditions(self, variables):
         c_s_xav = variables[
             "X-averaged " + self.domain.lower() + " particle concentration"
         ]
-        T_xav = pybamm.PrimaryBroadcast(
-            variables["X-averaged " + self.domain.lower() + " electrode temperature"],
-            c_s_xav.domain[0],
-        )
+        D_eff_xav = variables[
+            "X-averaged " + self.domain.lower() + " effective diffusivity"
+        ]
         j_xav = variables[
             "X-averaged "
             + self.domain.lower()
@@ -87,12 +89,7 @@ def set_boundary_conditions(self, variables):
         ]
 
         if self.domain == "Negative":
-            rbc = (
-                -self.param.C_n
-                * j_xav
-                / self.param.a_R_n
-                / pybamm.surf(self.param.D_n(c_s_xav, T_xav))
-            )
+            rbc = -self.param.C_n * j_xav / self.param.a_R_n / pybamm.surf(D_eff_xav)
 
         elif self.domain == "Positive":
             rbc = (
@@ -100,7 +97,7 @@ def set_boundary_conditions(self, variables):
                 * j_xav
                 / self.param.a_R_p
                 / self.param.gamma_p
-                / pybamm.surf(self.param.D_p(c_s_xav, T_xav))
+                / pybamm.surf(D_eff_xav)
             )
 
         self.boundary_conditions = {
diff --git a/pybamm/models/submodels/particle/no_distribution/x_averaged_polynomial_profile.py b/pybamm/models/submodels/particle/no_distribution/x_averaged_polynomial_profile.py
index c05ed3851e..356d3e797a 100644
--- a/pybamm/models/submodels/particle/no_distribution/x_averaged_polynomial_profile.py
+++ b/pybamm/models/submodels/particle/no_distribution/x_averaged_polynomial_profile.py
@@ -3,13 +3,13 @@
 #
 import pybamm
 
-from ..base_particle import BaseParticle
+from ..base_fickian import BaseFickian
 
 
-class XAveragedPolynomialProfile(BaseParticle):
+class XAveragedPolynomialProfile(BaseFickian):
     """
-    Class for molar conservation in a single x-averaged particle with
-    an assumed polynomial concentration profile in r. Model equations from [1]_.
+    Class for molar conservation in a single x-averaged particle employing Fick's law,
+    with an assumed polynomial concentration profile in r. Model equations from [1]_.
 
     Parameters
     ----------
@@ -20,6 +20,9 @@ class XAveragedPolynomialProfile(BaseParticle):
     name : str
         The name of the polynomial approximation to be used. Can be "uniform
         profile", "quadratic profile" or "quartic profile".
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     References
     ----------
@@ -30,8 +33,8 @@ class XAveragedPolynomialProfile(BaseParticle):
     **Extends:** :class:`pybamm.particle.BaseParticle`
     """
 
-    def __init__(self, param, domain, name):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, name, options):
+        super().__init__(param, domain, options)
         self.name = name
 
         pybamm.citations.register("Subramanian2005")
@@ -39,13 +42,13 @@ def __init__(self, param, domain, name):
     def get_fundamental_variables(self):
         # For all orders we solve an equation for the average concentration
         if self.domain == "Negative":
-            c_s_rxav = pybamm.standard_variables.c_s_n_rxav
+            c_s_av = pybamm.standard_variables.c_s_n_av
 
         elif self.domain == "Positive":
-            c_s_rxav = pybamm.standard_variables.c_s_p_rxav
+            c_s_av = pybamm.standard_variables.c_s_p_av
 
         variables = {
-            "Average " + self.domain.lower() + " particle concentration": c_s_rxav
+            "Average " + self.domain.lower() + " particle concentration": c_s_av
         }
 
         # For the fourth order polynomial approximation we also solve an
@@ -55,39 +58,37 @@ def get_fundamental_variables(self):
         # gradient q = dc/dr
         if self.name == "quartic profile":
             if self.domain == "Negative":
-                q_s_rxav = pybamm.standard_variables.q_s_n_rxav
+                q_s_av = pybamm.standard_variables.q_s_n_av
             elif self.domain == "Positive":
-                q_s_rxav = pybamm.standard_variables.q_s_p_rxav
+                q_s_av = pybamm.standard_variables.q_s_p_av
             variables.update(
                 {
                     "Average "
                     + self.domain.lower()
-                    + " particle concentration gradient": q_s_rxav
+                    + " particle concentration gradient": q_s_av
                 }
             )
 
         return variables
 
     def get_coupled_variables(self, variables):
-        c_s_rxav = variables[
-            "Average " + self.domain.lower() + " particle concentration"
-        ]
+        c_s_av = variables["Average " + self.domain.lower() + " particle concentration"]
+        T_av = (
+            variables["X-averaged " + self.domain.lower() + " electrode temperature"],
+        )
+        D_eff_av = self._get_effective_diffusivity(c_s_av, T_av)
         i_boundary_cc = variables["Current collector current density"]
         a_av = variables[
             "X-averaged "
             + self.domain.lower()
             + " electrode surface area to volume ratio"
         ]
-        T_xav = pybamm.PrimaryBroadcast(
-            variables["X-averaged " + self.domain.lower() + " electrode temperature"],
-            [self.domain.lower() + " particle"],
-        )
 
         # Set surface concentration based on polynomial order
         if self.name == "uniform profile":
             # The concentration is uniform so the surface value is equal to
             # the average
-            c_s_surf_xav = c_s_rxav
+            c_s_surf_xav = c_s_av
         elif self.name == "quadratic profile":
             # The surface concentration is computed from the average concentration
             # and boundary flux
@@ -106,72 +107,53 @@ def get_coupled_variables(self, variables):
             # than DAEs.
             if self.domain == "Negative":
                 j_xav = i_boundary_cc / (a_av * self.param.l_n)
-                c_s_surf_xav = c_s_rxav - self.param.C_n * (
-                    j_xav
-                    / 5
-                    / self.param.a_R_n
-                    / self.param.D_n(c_s_rxav, pybamm.surf(T_xav))
+                c_s_surf_xav = c_s_av - self.param.C_n * (
+                    j_xav / 5 / self.param.a_R_n / D_eff_av
                 )
 
             if self.domain == "Positive":
                 j_xav = -i_boundary_cc / (a_av * self.param.l_p)
-                c_s_surf_xav = c_s_rxav - self.param.C_p * (
-                    j_xav
-                    / 5
-                    / self.param.a_R_p
-                    / self.param.gamma_p
-                    / self.param.D_p(c_s_rxav, pybamm.surf(T_xav))
+                c_s_surf_xav = c_s_av - self.param.C_p * (
+                    j_xav / 5 / self.param.a_R_p / self.param.gamma_p / D_eff_av
                 )
         elif self.name == "quartic profile":
             # The surface concentration is computed from the average concentration,
             # the average concentration gradient and the boundary flux (see notes
-            # for the case order=2)
-            q_s_rxav = variables[
+            # for the quadratic profile)
+            q_s_av = variables[
                 "Average " + self.domain.lower() + " particle concentration gradient"
             ]
             if self.domain == "Negative":
                 j_xav = i_boundary_cc / (a_av * self.param.l_n)
                 c_s_surf_xav = (
-                    c_s_rxav
-                    + 8 * q_s_rxav / 35
-                    - self.param.C_n
-                    * (
-                        j_xav
-                        / 35
-                        / self.param.a_R_n
-                        / self.param.D_n(c_s_rxav, pybamm.surf(T_xav))
-                    )
+                    c_s_av
+                    + 8 * q_s_av / 35
+                    - self.param.C_n * (j_xav / 35 / self.param.a_R_n / D_eff_av)
                 )
 
             if self.domain == "Positive":
                 j_xav = -i_boundary_cc / (a_av * self.param.l_p)
                 c_s_surf_xav = (
-                    c_s_rxav
-                    + 8 * q_s_rxav / 35
+                    c_s_av
+                    + 8 * q_s_av / 35
                     - self.param.C_p
-                    * (
-                        j_xav
-                        / 35
-                        / self.param.a_R_p
-                        / self.param.gamma_p
-                        / self.param.D_p(c_s_rxav, pybamm.surf(T_xav))
-                    )
+                    * (j_xav / 35 / self.param.a_R_p / self.param.gamma_p / D_eff_av)
                 )
 
         # Set concentration depending on polynomial order
         if self.name == "uniform profile":
             # The concentration is uniform
             c_s_xav = pybamm.PrimaryBroadcast(
-                c_s_rxav, [self.domain.lower() + " particle"]
+                c_s_av, [self.domain.lower() + " particle"]
             )
         elif self.name == "quadratic profile":
             # The concentration is given by c = A + B*r**2
             A = pybamm.PrimaryBroadcast(
-                (1 / 2) * (5 * c_s_rxav - 3 * c_s_surf_xav),
+                (1 / 2) * (5 * c_s_av - 3 * c_s_surf_xav),
                 [self.domain.lower() + " particle"],
             )
             B = pybamm.PrimaryBroadcast(
-                (5 / 2) * (c_s_surf_xav - c_s_rxav), [self.domain.lower() + " particle"]
+                (5 / 2) * (c_s_surf_xav - c_s_av), [self.domain.lower() + " particle"]
             )
             if self.domain == "Negative":
                 # Since c_s_xav doesn't depend on x, we need to define a spatial
@@ -199,15 +181,15 @@ def get_coupled_variables(self, variables):
         elif self.name == "quartic profile":
             # The concentration is given by c = A + B*r**2 + C*r**4
             A = pybamm.PrimaryBroadcast(
-                39 * c_s_surf_xav / 4 - 3 * q_s_rxav - 35 * c_s_rxav / 4,
+                39 * c_s_surf_xav / 4 - 3 * q_s_av - 35 * c_s_av / 4,
                 [self.domain.lower() + " particle"],
             )
             B = pybamm.PrimaryBroadcast(
-                -35 * c_s_surf_xav + 10 * q_s_rxav + 35 * c_s_rxav,
+                -35 * c_s_surf_xav + 10 * q_s_av + 35 * c_s_av,
                 [self.domain.lower() + " particle"],
             )
             C = pybamm.PrimaryBroadcast(
-                105 * c_s_surf_xav / 4 - 7 * q_s_rxav - 105 * c_s_rxav / 4,
+                105 * c_s_surf_xav / 4 - 7 * q_s_av - 105 * c_s_av / 4,
                 [self.domain.lower() + " particle"],
             )
             if self.domain == "Negative":
@@ -239,6 +221,12 @@ def get_coupled_variables(self, variables):
         )
 
         # Set flux based on polynomial order
+        T_xav = pybamm.PrimaryBroadcast(
+            T_av,
+            [self.domain.lower() + " particle"],
+        )
+        D_eff_xav = self._get_effective_diffusivity(c_s_xav, T_xav)
+
         if self.name == "uniform profile":
             # The flux is zero since there is no concentration gradient
             N_s_xav = pybamm.FullBroadcastToEdges(
@@ -247,37 +235,36 @@ def get_coupled_variables(self, variables):
         elif self.name == "quadratic profile":
             # The flux may be computed directly from the polynomial for c
             if self.domain == "Negative":
-                N_s_xav = (
-                    -self.param.D_n(c_s_xav, T_xav) * 5 * (c_s_surf_xav - c_s_rxav) * r
-                )
+                N_s_xav = -D_eff_xav * 5 * (c_s_surf_xav - c_s_av) * r
             if self.domain == "Positive":
-                N_s_xav = (
-                    -self.param.D_p(c_s_xav, T_xav) * 5 * (c_s_surf_xav - c_s_rxav) * r
-                )
+                N_s_xav = -D_eff_xav * 5 * (c_s_surf_xav - c_s_av) * r
         elif self.name == "quartic profile":
-            q_s_rxav = variables[
+            q_s_av = variables[
                 "Average " + self.domain.lower() + " particle concentration gradient"
             ]
             # The flux may be computed directly from the polynomial for c
             if self.domain == "Negative":
-                N_s_xav = -self.param.D_n(c_s_xav, T_xav) * (
-                    (-70 * c_s_surf_xav + 20 * q_s_rxav + 70 * c_s_rxav) * r
-                    + (105 * c_s_surf_xav - 28 * q_s_rxav - 105 * c_s_rxav) * r ** 3
+                N_s_xav = -D_eff_xav * (
+                    (-70 * c_s_surf_xav + 20 * q_s_av + 70 * c_s_av) * r
+                    + (105 * c_s_surf_xav - 28 * q_s_av - 105 * c_s_av) * r ** 3
                 )
             elif self.domain == "Positive":
-                N_s_xav = -self.param.D_p(c_s_xav, T_xav) * (
-                    (-70 * c_s_surf_xav + 20 * q_s_rxav + 70 * c_s_rxav) * r
-                    + (105 * c_s_surf_xav - 28 * q_s_rxav - 105 * c_s_rxav) * r ** 3
+                N_s_xav = -D_eff_xav * (
+                    (-70 * c_s_surf_xav + 20 * q_s_av + 70 * c_s_av) * r
+                    + (105 * c_s_surf_xav - 28 * q_s_av - 105 * c_s_av) * r ** 3
                 )
 
+        D_eff = pybamm.SecondaryBroadcast(
+            D_eff_xav, [self._domain.lower() + " electrode"]
+        )
         N_s = pybamm.SecondaryBroadcast(N_s_xav, [self._domain.lower() + " electrode"])
 
         variables.update(
             self._get_standard_concentration_variables(
-                c_s, c_s_av=c_s_rxav, c_s_surf=c_s_surf
+                c_s, c_s_av=c_s_av, c_s_surf=c_s_surf
             )
         )
-        variables.update(self._get_standard_flux_variables(N_s, N_s_xav))
+        variables.update(self._get_standard_flux_variables(N_s, N_s_xav, D_eff))
         variables.update(self._get_total_concentration_variables(variables))
 
         return variables
@@ -288,9 +275,7 @@ def set_rhs(self, variables):
         # using this model with 2D current collectors with the finite element
         # method (see #1399)
 
-        c_s_rxav = variables[
-            "Average " + self.domain.lower() + " particle concentration"
-        ]
+        c_s_av = variables["Average " + self.domain.lower() + " particle concentration"]
         j_xav = variables[
             "X-averaged "
             + self.domain.lower()
@@ -298,51 +283,41 @@ def set_rhs(self, variables):
         ]
 
         if self.domain == "Negative":
-            self.rhs = {
-                c_s_rxav: pybamm.source(-3 * j_xav / self.param.a_R_n, c_s_rxav)
-            }
+            self.rhs = {c_s_av: pybamm.source(-3 * j_xav / self.param.a_R_n, c_s_av)}
 
         elif self.domain == "Positive":
             self.rhs = {
-                c_s_rxav: pybamm.source(
-                    -3 * j_xav / self.param.a_R_p / self.param.gamma_p, c_s_rxav
+                c_s_av: pybamm.source(
+                    -3 * j_xav / self.param.a_R_p / self.param.gamma_p, c_s_av
                 )
             }
 
         if self.name == "quartic profile":
             # We solve an extra ODE for the average particle concentration gradient
-            q_s_rxav = variables[
+            q_s_av = variables[
                 "Average " + self.domain.lower() + " particle concentration gradient"
             ]
-            c_s_surf_xav = variables[
-                "X-averaged " + self.domain.lower() + " particle surface concentration"
-            ]
-            T_xav = variables[
-                "X-averaged " + self.domain.lower() + " electrode temperature"
+            D_eff_xav = variables[
+                "X-averaged " + self.domain.lower() + " effective diffusivity"
             ]
+
             if self.domain == "Negative":
                 self.rhs.update(
                     {
-                        q_s_rxav: pybamm.source(
-                            -30
-                            * self.param.D_n(c_s_surf_xav, T_xav)
-                            * q_s_rxav
-                            / self.param.C_n
+                        q_s_av: pybamm.source(
+                            -30 * pybamm.surf(D_eff_xav) * q_s_av / self.param.C_n
                             - 45 * j_xav / self.param.a_R_n / 2,
-                            q_s_rxav,
+                            q_s_av,
                         )
                     }
                 )
             elif self.domain == "Positive":
                 self.rhs.update(
                     {
-                        q_s_rxav: pybamm.source(
-                            -30
-                            * self.param.D_p(c_s_surf_xav, T_xav)
-                            * q_s_rxav
-                            / self.param.C_p
+                        q_s_av: pybamm.source(
+                            -30 * pybamm.surf(D_eff_xav) * q_s_av / self.param.C_p
                             - 45 * j_xav / self.param.a_R_p / self.param.gamma_p / 2,
-                            q_s_rxav,
+                            q_s_av,
                         )
                     }
                 )
@@ -353,9 +328,7 @@ def set_initial_conditions(self, variables):
         so we arbitrarily evaluate them at x=0 in the negative electrode and x=1 in the
         positive electrode (they will usually be constant)
         """
-        c_s_rxav = variables[
-            "Average " + self.domain.lower() + " particle concentration"
-        ]
+        c_s_av = variables["Average " + self.domain.lower() + " particle concentration"]
 
         if self.domain == "Negative":
             c_init = self.param.c_n_init(0)
@@ -363,11 +336,11 @@ def set_initial_conditions(self, variables):
         elif self.domain == "Positive":
             c_init = self.param.c_p_init(1)
 
-        self.initial_conditions = {c_s_rxav: c_init}
+        self.initial_conditions = {c_s_av: c_init}
         if self.name == "quartic profile":
             # We also need to provide an initial condition for the average
             # concentration gradient
-            q_s_rxav = variables[
+            q_s_av = variables[
                 "Average " + self.domain.lower() + " particle concentration gradient"
             ]
-            self.initial_conditions.update({q_s_rxav: 0})
+            self.initial_conditions.update({q_s_av: 0})
diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py
index 4bc8f16c4b..82540327fe 100644
--- a/pybamm/parameters/lithium_ion_parameters.py
+++ b/pybamm/parameters/lithium_ion_parameters.py
@@ -245,27 +245,27 @@ def _set_dimensional_parameters(self):
             "Initial concentration in electrolyte [mol.m-3]"
         )
 
-        # mechanical parameters
-        self.nu_p = pybamm.Parameter("Positive electrode Poisson's ratio")
-        self.E_p = pybamm.Parameter("Positive electrode Young's modulus [Pa]")
-        self.c_p_0_dim = pybamm.Parameter(
-            "Positive electrode reference concentration for free of deformation "
-            "[mol.m-3]"
-        )
-        self.Omega_p = pybamm.Parameter(
-            "Positive electrode partial molar volume [m3.mol-1]"
-        )
+        # Mechanical parameters
         self.nu_n = pybamm.Parameter("Negative electrode Poisson's ratio")
+        self.nu_p = pybamm.Parameter("Positive electrode Poisson's ratio")
         self.E_n = pybamm.Parameter("Negative electrode Young's modulus [Pa]")
+        self.E_p = pybamm.Parameter("Positive electrode Young's modulus [Pa]")
         self.c_n_0_dim = pybamm.Parameter(
             "Negative electrode reference concentration for free of deformation "
             "[mol.m-3]"
         )
+        self.c_p_0_dim = pybamm.Parameter(
+            "Positive electrode reference concentration for free of deformation "
+            "[mol.m-3]"
+        )
         self.Omega_n = pybamm.Parameter(
             "Negative electrode partial molar volume [m3.mol-1]"
         )
-        self.l_cr_p_0 = pybamm.Parameter("Positive electrode initial crack length [m]")
+        self.Omega_p = pybamm.Parameter(
+            "Positive electrode partial molar volume [m3.mol-1]"
+        )
         self.l_cr_n_0 = pybamm.Parameter("Negative electrode initial crack length [m]")
+        self.l_cr_p_0 = pybamm.Parameter("Positive electrode initial crack length [m]")
         self.w_cr = pybamm.Parameter("Negative electrode initial crack width [m]")
         self.rho_cr_n_dim = pybamm.Parameter(
             "Negative electrode number of cracks per unit area [m-2]"
@@ -274,27 +274,25 @@ def _set_dimensional_parameters(self):
             "Positive electrode number of cracks per unit area [m-2]"
         )
         self.b_cr_n = pybamm.Parameter("Negative electrode Paris' law constant b")
+        self.b_cr_p = pybamm.Parameter("Positive electrode Paris' law constant b")
         self.m_cr_n = pybamm.Parameter("Negative electrode Paris' law constant m")
+        self.m_cr_p = pybamm.Parameter("Positive electrode Paris' law constant m")
         self.Eac_cr_n = pybamm.Parameter(
             "Negative electrode activation energy for cracking rate [kJ.mol-1]"
         )  # noqa
-        self.b_cr_p = pybamm.Parameter("Positive electrode Paris' law constant b")
-        self.m_cr_p = pybamm.Parameter("Positive electrode Paris' law constant m")
         self.Eac_cr_p = pybamm.Parameter(
             "Positive electrode activation energy for cracking rate [kJ.mol-1]"
         )  # noqa
-        self.alpha_T_cell_dim = pybamm.Parameter(
-            "Cell thermal expansion coefficient [m.K-1]"
+        # intermediate variables  [K*m^3/mol]
+        self.theta_n_dim = (
+            (self.Omega_n / self.R) * 2 * self.Omega_n * self.E_n / 9 / (1 - self.nu_n)
         )
-        self.R_const = pybamm.constants.R
         self.theta_p_dim = (
-            self.Omega_p ** 2 / self.R_const * 2 / 9 * self.E_p / (1 - self.nu_p)
+            (self.Omega_p / self.R) * 2 * self.Omega_p * self.E_p / 9 / (1 - self.nu_p)
         )
-        # intermediate variable  [K*m^3/mol]
-        self.theta_n_dim = (
-            self.Omega_n ** 2 / self.R_const * 2 / 9 * self.E_n / (1 - self.nu_n)
+        self.alpha_T_cell_dim = pybamm.Parameter(
+            "Cell thermal expansion coefficient [m.K-1]"
         )
-        # intermediate variable  [K*m^3/mol]
 
         # Electrode capacities
         x_n = pybamm.SpatialVariable(
@@ -330,7 +328,8 @@ def _set_dimensional_parameters(self):
 
         self.n_Li_particles_init = self.n_Li_n_init + self.n_Li_p_init
         self.n_Li_init = self.n_Li_particles_init + self.n_Li_e_init
-        # loss of active material parameters
+
+        # Loss of active material parameters
         self.m_LAM_n = pybamm.Parameter(
             "Negative electrode LAM constant exponential term"
         )
@@ -384,32 +383,16 @@ def D_n_dimensional(self, sto, T):
         """Dimensional diffusivity in negative particle. Note this is defined as a
         function of stochiometry"""
         inputs = {"Negative particle stoichiometry": sto, "Temperature [K]": T}
-        crack = self.options["particle mechanics"]
-        if crack != "none" or (isinstance(crack, tuple) and crack[0] != "none"):
-            mech_effects = (
-                1 + self.theta_n_dim * (sto * self.c_n_max - self.c_n_0_dim) / T
-            )
-        else:
-            mech_effects = 1
-        return (
-            pybamm.FunctionParameter("Negative electrode diffusivity [m2.s-1]", inputs)
-            * mech_effects
+        return pybamm.FunctionParameter(
+            "Negative electrode diffusivity [m2.s-1]", inputs
         )
 
     def D_p_dimensional(self, sto, T):
         """Dimensional diffusivity in positive particle. Note this is defined as a
         function of stochiometry"""
         inputs = {"Positive particle stoichiometry": sto, "Temperature [K]": T}
-        crack = self.options["particle mechanics"]
-        if crack != "none" or (isinstance(crack, tuple) and crack[1] != "none"):
-            mech_effects = (
-                1 + self.theta_p_dim * (sto * self.c_p_max - self.c_p_0_dim) / T
-            )
-        else:
-            mech_effects = 1
-        return (
-            pybamm.FunctionParameter("Positive electrode diffusivity [m2.s-1]", inputs)
-            * mech_effects
+        return pybamm.FunctionParameter(
+            "Positive electrode diffusivity [m2.s-1]", inputs
         )
 
     def j0_n_dimensional(self, c_e, c_s_surf, T):
@@ -854,8 +837,8 @@ def _set_dimensionless_parameters(self):
         # Dimensionless mechanical parameters
         self.rho_cr_n = self.rho_cr_n_dim * self.l_cr_n_0 * self.w_cr
         self.rho_cr_p = self.rho_cr_p_dim * self.l_cr_p_0 * self.w_cr
-        self.theta_p = self.theta_p_dim * self.c_p_max / self.Delta_T
-        self.theta_n = self.theta_n_dim * self.c_n_max / self.Delta_T
+        self.theta_p = self.theta_p_dim * self.c_p_max / self.T_ref
+        self.theta_n = self.theta_n_dim * self.c_n_max / self.T_ref
         self.c_p_0 = self.c_p_0_dim / self.c_p_max
         self.c_n_0 = self.c_n_0_dim / self.c_n_max
         self.t0_cr = 3600 / self.C_rate / self.timescale
@@ -1040,25 +1023,9 @@ def rho(self, T):
             + self.rho_cp(T) * self.l_cp
         ) / self.l
 
-    def _set_input_current(self):
-        """Set the input current"""
-
-        self.dimensional_current_with_time = pybamm.FunctionParameter(
-            "Current function [A]", {"Time [s]": pybamm.t * self.timescale}
-        )
-        self.dimensional_current_density_with_time = (
-            self.dimensional_current_with_time
-            / (self.n_electrodes_parallel * self.geo.A_cc)
-        )
-        self.current_with_time = (
-            self.dimensional_current_with_time
-            / self.I_typ
-            * pybamm.Function(np.sign, self.I_typ)
-        )
-
     def t_n_change(self, sto):
         """
-        Dimentionless volume change for the negative electrode;
+        Dimensionless volume change for the negative electrode;
         sto should be R-averaged
         """
         return pybamm.FunctionParameter(
@@ -1067,7 +1034,7 @@ def t_n_change(self, sto):
 
     def t_p_change(self, sto):
         """
-        Dimentionless volume change for the positive electrode;
+        Dimensionless volume change for the positive electrode;
         sto should be R-averaged
         """
         return pybamm.FunctionParameter(
@@ -1076,7 +1043,7 @@ def t_p_change(self, sto):
 
     def k_cr_p(self, T):
         """
-        Dimentionless cracking rate for the positive electrode;
+        Dimensionless cracking rate for the positive electrode;
         """
         T_dim = self.Delta_T * T + self.T_ref
         delta_k_cr = self.E_p ** self.m_cr_p * self.l_cr_p_0 ** (self.m_cr_p / 2 - 1)
@@ -1089,7 +1056,7 @@ def k_cr_p(self, T):
 
     def k_cr_n(self, T):
         """
-        Dimentionless cracking rate for the negative electrode;
+        Dimensionless cracking rate for the negative electrode;
         """
         T_dim = self.Delta_T * T + self.T_ref
         delta_k_cr = self.E_n ** self.m_cr_n * self.l_cr_n_0 ** (self.m_cr_n / 2 - 1)
@@ -1100,6 +1067,22 @@ def k_cr_n(self, T):
             * delta_k_cr
         )
 
+    def _set_input_current(self):
+        """Set the input current"""
+
+        self.dimensional_current_with_time = pybamm.FunctionParameter(
+            "Current function [A]", {"Time [s]": pybamm.t * self.timescale}
+        )
+        self.dimensional_current_density_with_time = (
+            self.dimensional_current_with_time
+            / (self.n_electrodes_parallel * self.geo.A_cc)
+        )
+        self.current_with_time = (
+            self.dimensional_current_with_time
+            / self.I_typ
+            * pybamm.Function(np.sign, self.I_typ)
+        )
+
     @property
     def options(self):
         return self._options