diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8f8ff3599a..ed5713e5c6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -87,6 +87,8 @@
 - Fix to simulate c_rate steps with drive cycles ([#3186](https://github.com/pybamm-team/PyBaMM/pull/3186))
 - Always save last cycle in experiment, to fix issues with `starting_solution` and `last_state` ([#3177](https://github.com/pybamm-team/PyBaMM/pull/3177))
 - Fix simulations with `starting_solution` to work with `start_time` experiments ([#3177](https://github.com/pybamm-team/PyBaMM/pull/3177))
+- Fixed how `Loss of lithium due to loss of active material` variable is calculated ([#3171](https://github.com/pybamm-team/PyBaMM/pull/3171))
+- Rewritten SEI submodel so that concentration, not thickness, is the fundamental variable. This change ensures lithium is always conserved. ([#3171](https://github.com/pybamm-team/PyBaMM/pull/3171))
 - Fix SEI Example Notebook ([#3166](https://github.com/pybamm-team/PyBaMM/pull/3166))
 - Thevenin() model is now constructed with standard variables: `Time [s]`, `Time [min]`, `Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
 - Error generated when invalid parameter values are passed ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
diff --git a/docs/source/examples/notebooks/models/coupled-degradation.ipynb b/docs/source/examples/notebooks/models/coupled-degradation.ipynb
index 8c083d986a..901a862b1d 100644
--- a/docs/source/examples/notebooks/models/coupled-degradation.ipynb
+++ b/docs/source/examples/notebooks/models/coupled-degradation.ipynb
@@ -57,7 +57,6 @@
     "        \"particle mechanics\": (\"swelling and cracking\", \"swelling only\"),\n",
     "        \"SEI on cracks\": \"true\",\n",
     "        \"loss of active material\": \"stress-driven\",\n",
-    "        \"calculate discharge energy\": \"true\",  # for compatibility with older PyBaMM versions\n",
     "    }\n",
     ")"
    ]
diff --git a/docs/source/examples/notebooks/models/using-submodels.ipynb b/docs/source/examples/notebooks/models/using-submodels.ipynb
index 5c00de05da..97ddc563be 100644
--- a/docs/source/examples/notebooks/models/using-submodels.ipynb
+++ b/docs/source/examples/notebooks/models/using-submodels.ipynb
@@ -515,15 +515,27 @@
     "model.submodels[\"Negative sei\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"negative\", model.options\n",
     ")\n",
+    "model.submodels[\"Negative sei thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"negative\", model.options\n",
+    ")\n",
     "model.submodels[\"Positive sei\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"positive\", model.options\n",
     ")\n",
+    "model.submodels[\"Positive sei thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"positive\", model.options\n",
+    ")\n",
     "model.submodels[\"Negative sei on cracks\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"negative\", model.options, cracks=True\n",
     ")\n",
+    "model.submodels[\"Negative sei on cracks thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"negative\", model.options, cracks=True\n",
+    ")\n",
     "model.submodels[\"Positive sei on cracks\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"positive\", model.options, cracks=True\n",
     ")\n",
+    "model.submodels[\"Positive sei on cracks thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"positive\", model.options, cracks=True\n",
+    ")\n",
     "model.submodels[\"Negative lithium plating\"] = pybamm.lithium_plating.NoPlating(\n",
     "    model.param, \"Negative\"\n",
     ")\n",
diff --git a/examples/scripts/custom_model.py b/examples/scripts/custom_model.py
index 4c6dcf3fad..6506137470 100644
--- a/examples/scripts/custom_model.py
+++ b/examples/scripts/custom_model.py
@@ -75,9 +75,15 @@
     model.submodels[f"{domain} sei"] = pybamm.sei.NoSEI(
         model.param, domain, model.options
     )
+    model.submodels[f"{domain} sei thickness"] = pybamm.sei.SEIThickness(
+        model.param, domain, model.options
+    )
     model.submodels[f"{domain} sei on cracks"] = pybamm.sei.NoSEI(
         model.param, domain, model.options, cracks=True
     )
+    model.submodels[f"{domain} sei on cracks thickness"] = pybamm.sei.SEIThickness(
+        model.param, domain, model.options, cracks=True
+    )
     model.submodels[f"{domain} lithium plating"] = pybamm.lithium_plating.NoPlating(
         model.param, domain
     )
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 fbe19b0d42..6262acf22e 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
@@ -256,12 +256,6 @@ def set_open_circuit_potential_submodel(self):
 
     def set_sei_submodel(self):
         for domain in ["negative", "positive"]:
-            if self.options.electrode_types[domain] == "planar":
-                reaction_loc = "interface"
-            elif self.options["x-average side reactions"] == "true":
-                reaction_loc = "x-average"
-            else:
-                reaction_loc = "full electrode"
             sei_option = getattr(self.options, domain)["SEI"]
             phases = self.options.phases[domain]
             for phase in phases:
@@ -275,12 +269,16 @@ def set_sei_submodel(self):
                     submodel = pybamm.sei.SEIGrowth(
                         self.param,
                         domain,
-                        reaction_loc,
                         self.options,
                         phase,
                         cracks=False,
                     )
                 self.submodels[f"{domain} {phase} sei"] = submodel
+                self.submodels[
+                    f"{domain} {phase} sei thickness"
+                ] = pybamm.sei.SEIThickness(
+                    self.param, domain, self.options, phase, cracks=False
+                )
             if len(phases) > 1:
                 self.submodels[f"{domain} total sei"] = pybamm.sei.TotalSEI(
                     self.param, domain, self.options
@@ -304,19 +302,23 @@ def set_sei_on_cracks_submodel(self):
                             self.param, domain, self.options, phase, cracks=True
                         )
                     else:
-                        if self.options["x-average side reactions"] == "true":
-                            reaction_loc = "x-average"
-                        else:
-                            reaction_loc = "full electrode"
                         submodel = pybamm.sei.SEIGrowth(
                             self.param,
                             domain,
-                            reaction_loc,
                             self.options,
                             phase,
                             cracks=True,
                         )
                     self.submodels[f"{domain} {phase} sei on cracks"] = submodel
+                    self.submodels[
+                        f"{domain} {phase} sei on cracks thickness"
+                    ] = pybamm.sei.SEIThickness(
+                        self.param,
+                        domain,
+                        self.options,
+                        phase,
+                        cracks=True,
+                    )
                 if len(phases) > 1:
                     self.submodels[
                         f"{domain} total sei on cracks"
diff --git a/pybamm/models/submodels/active_material/loss_active_material.py b/pybamm/models/submodels/active_material/loss_active_material.py
index 7816122e07..f21423194b 100644
--- a/pybamm/models/submodels/active_material/loss_active_material.py
+++ b/pybamm/models/submodels/active_material/loss_active_material.py
@@ -153,12 +153,12 @@ def set_rhs(self, variables):
             f"in {domain} electrode [mol]"
         ]
         # Multiply by mol.m-3 * m3 to get mol
-        c_s_av = variables[f"Average {domain} particle concentration [mol.m-3]"]
+        c_s_rav = variables[f"R-averaged {domain} particle concentration [mol.m-3]"]
         V = self.domain_param.L * self.param.A_cc
 
         self.rhs = {
             # minus sign because eps_solid is decreasing and LLI measures positive
-            lli_due_to_lam: -c_s_av * V * pybamm.x_average(deps_solid_dt),
+            lli_due_to_lam: -V * pybamm.x_average(c_s_rav * deps_solid_dt),
             eps_solid: deps_solid_dt,
         }
 
diff --git a/pybamm/models/submodels/interface/sei/__init__.py b/pybamm/models/submodels/interface/sei/__init__.py
index 5f151daf0e..a2166ba6d1 100644
--- a/pybamm/models/submodels/interface/sei/__init__.py
+++ b/pybamm/models/submodels/interface/sei/__init__.py
@@ -3,3 +3,4 @@
 from .no_sei import NoSEI
 from .constant_sei import ConstantSEI
 from .sei_growth import SEIGrowth
+from .sei_thickness import SEIThickness
diff --git a/pybamm/models/submodels/interface/sei/base_sei.py b/pybamm/models/submodels/interface/sei/base_sei.py
index b0e8db56c6..168556c056 100644
--- a/pybamm/models/submodels/interface/sei/base_sei.py
+++ b/pybamm/models/submodels/interface/sei/base_sei.py
@@ -55,17 +55,17 @@ def get_coupled_variables(self, variables):
 
         return variables
 
-    def _get_standard_thickness_variables(self, L_inner, L_outer):
+    def _get_standard_concentration_variables(self, c_inner, c_outer):
         """
         A private function to obtain the standard variables which
-        can be derived from the local SEI thickness.
+        can be derived from the local SEI concentration.
 
         Parameters
         ----------
-        L_inner : :class:`pybamm.Symbol`
-            The inner SEI thickness.
-        L_outer : :class:`pybamm.Symbol`
-            The outer SEI thickness.
+        c_inner : :class:`pybamm.Symbol`
+            The inner SEI concentration.
+        c_outer : :class:`pybamm.Symbol`
+            The outer SEI concentration.
 
         Returns
         -------
@@ -73,199 +73,47 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
             The variables which can be derived from the SEI thicknesses.
         """
         domain, Domain = self.domain_Domain
-        variables = {
-            f"{Domain} inner {self.reaction_name}thickness [m]": L_inner,
-            f"{Domain} outer {self.reaction_name}thickness [m]": L_outer,
-        }
-
-        if self.reaction_loc != "interface":
-            L_inner_av = pybamm.x_average(L_inner)
-            L_outer_av = pybamm.x_average(L_outer)
-            variables.update(
-                {
-                    f"X-averaged {domain} inner {self.reaction_name}"
-                    "thickness [m]": L_inner_av,
-                    f"X-averaged {domain} outer {self.reaction_name}"
-                    "thickness [m]": L_outer_av,
-                }
-            )
-        # Get variables related to the total thickness
-        L_sei = L_inner + L_outer
-        variables.update(self._get_standard_total_thickness_variables(L_sei))
-
-        return variables
-
-    def _get_standard_total_thickness_variables(self, L_sei):
-        """Update variables related to total SEI thickness."""
-        domain, Domain = self.domain_Domain
+        if self.reaction_loc == "interface":  # c is an interfacial quantity [mol.m-2]
+            variables = {
+                f"{Domain} inner {self.reaction_name}concentration [mol.m-2]": c_inner,
+                f"{Domain} outer {self.reaction_name}concentration [mol.m-2]": c_outer,
+            }
+        else:  # c is a bulk quantity [mol.m-3]
+            c_inner_av = pybamm.x_average(c_inner)
+            c_outer_av = pybamm.x_average(c_outer)
+            variables = {
+                f"{Domain} inner {self.reaction_name}concentration [mol.m-3]": c_inner,
+                f"{Domain} outer {self.reaction_name}concentration [mol.m-3]": c_outer,
+                f"X-averaged {domain} inner {self.reaction_name}"
+                "concentration [mol.m-3]": c_inner_av,
+                f"X-averaged {domain} outer {self.reaction_name}"
+                "concentration [mol.m-3]": c_outer_av,
+            }
+
+        # Get variables related to the total concentration
+        c_sei = c_inner + c_outer
+        variables.update(self._get_standard_total_concentration_variables(c_sei))
 
-        if isinstance(self, pybamm.sei.NoSEI):
-            R_sei = 1
-        else:
-            R_sei = self.phase_param.R_sei
-
-        variables = {
-            f"{Domain} {self.reaction_name}[m]": L_sei,
-            f"{Domain} total {self.reaction_name}thickness [m]": L_sei,
-        }
-        if self.reaction_loc != "interface":
-            L_sei_av = pybamm.x_average(L_sei)
-            variables.update(
-                {
-                    f"X-averaged {domain} {self.reaction_name}thickness [m]": L_sei_av,
-                    f"X-averaged {domain} total {self.reaction_name}"
-                    "thickness [m]": L_sei_av,
-                }
-            )
-            if self.reaction == "SEI":
-                variables.update(
-                    {
-                        f"X-averaged {domain} electrode resistance "
-                        "[Ohm.m2]": L_sei_av * R_sei,
-                    }
-                )
         return variables
 
-    def _get_standard_concentration_variables(self, variables):
-        """Update variables related to the SEI concentration."""
+    def _get_standard_total_concentration_variables(self, c_sei):
+        """Update variables related to total SEI concentration."""
         domain, Domain = self.domain_Domain
-        phase_param = self.phase_param
-        reaction_name = self.reaction_name
-
-        # 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):
-            L_to_n_inner = 0
-            L_to_n_outer = 0
-            n_SEI_0 = 0
-            n_crack_0 = 0
-            z_sei = 1
-        else:
-            if self.reaction_loc == "interface":
-                # m * (mol/m3) = mol/m2 (n is an interfacial quantity)
-                L_to_n_inner = 1 / phase_param.V_bar_inner
-                L_to_n_outer = 1 / phase_param.V_bar_outer
-                L_to_n_inner_0 = L_to_n_inner
-                L_to_n_outer_0 = L_to_n_outer
-            else:
-                # m * (mol/m4) = mol/m3 (n is a bulk quantity)
-                a = variables[
-                    f"{Domain} electrode {self.phase_name}"
-                    "surface area to volume ratio [m-1]"
-                ]
-                L_to_n_inner = a / phase_param.V_bar_inner
-                L_to_n_outer = a / phase_param.V_bar_outer
-                L_to_n_inner_0 = phase_param.a_typ / phase_param.V_bar_inner
-                L_to_n_outer_0 = phase_param.a_typ / phase_param.V_bar_outer
-            z_sei = phase_param.z_sei
-            L_inner_0 = phase_param.L_inner_0
-            L_outer_0 = phase_param.L_outer_0
-            L_inner_crack_0 = phase_param.L_inner_crack_0
-            L_outer_crack_0 = phase_param.L_outer_crack_0
-            n_SEI_0 = L_inner_0 * L_to_n_inner_0 + L_outer_0 * L_to_n_outer_0
-            n_crack_0 = (
-                L_inner_crack_0 * L_to_n_inner_0 + L_outer_crack_0 * L_to_n_outer_0
-            )
-
-        if self.reaction == "SEI":
-            L_inner = variables[f"{Domain} inner {reaction_name}thickness [m]"]
-            L_outer = variables[f"{Domain} outer {reaction_name}thickness [m]"]
-
-            n_inner = L_inner * L_to_n_inner  # inner SEI concentration
-            n_outer = L_outer * L_to_n_outer  # outer SEI concentration
-
-            n_inner_av = pybamm.x_average(n_inner)
-            n_outer_av = pybamm.x_average(n_outer)
-
-            n_SEI = n_inner + n_outer  # SEI concentration
-            n_SEI_xav = pybamm.x_average(n_SEI)
-            n_SEI_av = pybamm.yz_average(n_SEI_xav)
-
-            # Calculate change in SEI concentration with respect to initial state
-            delta_n_SEI = n_SEI_av - n_SEI_0
-
-            # Q_sei in mol
-            if self.reaction_loc == "interface":
-                L_k = 1
-            elif domain == "negative":
-                L_k = self.param.n.L
-            elif domain == "positive":
-                L_k = self.param.p.L
-
-            # Multiply delta_n_SEI by V_k to get total moles of SEI formed
-            # multiply by z_sei to get total lithium moles consumed by SEI
-            V_k = L_k * self.param.L_y * self.param.L_z
-            Q_sei = z_sei * delta_n_SEI * V_k
-
-            variables.update(
-                {
-                    f"{Domain} inner {reaction_name}concentration [mol.m-3]": n_inner,
-                    f"X-averaged {domain} inner {reaction_name}"
-                    "concentration [mol.m-3]": n_inner_av,
-                    f"{Domain} outer {reaction_name}concentration [mol.m-3]": n_outer,
-                    f"X-averaged {domain} outer {reaction_name}"
-                    "concentration [mol.m-3]": n_outer_av,
-                    f"{Domain} {reaction_name}concentration [mol.m-3]": n_SEI,
-                    f"X-averaged {domain} {reaction_name}"
-                    "concentration [mol.m-3]": n_SEI_xav,
-                    f"Loss of lithium to {domain} {reaction_name}[mol]": Q_sei,
-                    f"Loss of capacity to {domain} {reaction_name}[A.h]": Q_sei
-                    * self.param.F
-                    / 3600,
-                }
-            )
-        # Concentration variables are handled slightly differently for SEI on cracks
-        elif self.reaction == "SEI on cracks":
-            L_inner_cr = variables[f"{Domain} inner {reaction_name}thickness [m]"]
-            L_outer_cr = variables[f"{Domain} outer {reaction_name}thickness [m]"]
-            roughness = variables[f"{Domain} electrode roughness ratio"]
-
-            n_inner_cr = L_inner_cr * L_to_n_inner * (roughness - 1)
-            n_outer_cr = L_outer_cr * L_to_n_outer * (roughness - 1)
-
-            n_inner_cr_av = pybamm.x_average(n_inner_cr)
-            n_outer_cr_av = pybamm.x_average(n_outer_cr)
-
-            n_SEI_cr = n_inner_cr + n_outer_cr  # SEI on cracks concentration
-            n_SEI_cr_xav = pybamm.x_average(n_SEI_cr)
-            n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav)
-
-            # Calculate change in SEI cracks concentration
-            # Initial state depends on roughness (to avoid division by zero)
-            roughness_av = pybamm.yz_average(pybamm.x_average(roughness))
-            # choose an initial condition that is as close to zero to get the
-            # physics right, but doesn't cause a division by zero error
-            n_SEI_cr_init = n_crack_0 * (roughness_av - 1)
-            delta_n_SEI_cr = n_SEI_cr_av - n_SEI_cr_init
-
-            if domain == "negative":
-                L_k = self.param.n.L
-            elif domain == "positive":
-                L_k = self.param.p.L
-
-            # Q_sei_cr in mol
-            Q_sei_cr = z_sei * delta_n_SEI_cr * L_k * self.param.L_y * self.param.L_z
-
-            variables.update(
-                {
-                    f"{Domain} inner {reaction_name}"
-                    "concentration [mol.m-3]": n_inner_cr,
-                    f"X-averaged {domain} inner {reaction_name}"
-                    "concentration [mol.m-3]": n_inner_cr_av,
-                    f"{Domain} outer {reaction_name}"
-                    "concentration [mol.m-3]": n_outer_cr,
-                    f"X-averaged {domain} outer {reaction_name}"
-                    "concentration [mol.m-3]": n_outer_cr_av,
-                    f"{Domain} {reaction_name}" "concentration [mol.m-3]": n_SEI_cr,
-                    f"X-averaged {domain} {reaction_name}"
-                    "concentration [mol.m-3]": n_SEI_cr_xav,
-                    f"Loss of lithium to {domain} {reaction_name}[mol]": Q_sei_cr,
-                    f"Loss of capacity to {domain} {reaction_name}[A.h]": Q_sei_cr
-                    * self.param.F
-                    / 3600,
-                }
-            )
-
+        if self.reaction_loc == "interface":  # c is an interfacial quantity [mol.m-2]
+            variables = {
+                f"{Domain} {self.reaction_name}concentration [mol.m-2]": c_sei,
+                f"{Domain} total {self.reaction_name}concentration [mol.m-2]": c_sei,
+            }
+        else:  # c is a bulk quantity [mol.m-3]
+            c_xav = pybamm.x_average(c_sei)
+            variables = {
+                f"{Domain} {self.reaction_name}concentration [mol.m-3]": c_sei,
+                f"{Domain} total {self.reaction_name}concentration [mol.m-3]": c_sei,
+                f"X-averaged {domain} {self.reaction_name}"
+                "concentration [mol.m-3]": c_xav,
+                f"X-averaged {domain} total {self.reaction_name}"
+                "concentration [mol.m-3]": c_xav,
+            }
         return variables
 
     def _get_standard_reaction_variables(self, j_inner, j_outer):
diff --git a/pybamm/models/submodels/interface/sei/constant_sei.py b/pybamm/models/submodels/interface/sei/constant_sei.py
index 4c507eec3a..06fc9b1098 100644
--- a/pybamm/models/submodels/interface/sei/constant_sei.py
+++ b/pybamm/models/submodels/interface/sei/constant_sei.py
@@ -7,7 +7,7 @@
 
 class ConstantSEI(BaseModel):
     """
-    Class for SEI with constant thickness.
+    Class for SEI with constant concentration.
 
     Note that there is no SEI current, so we don't need to update the "sum of
     interfacial current densities" variables from
@@ -31,11 +31,19 @@ def __init__(self, param, domain, options, phase="primary"):
             self.reaction_loc = "full electrode"
 
     def get_fundamental_variables(self):
+        # Constant concentrations
         domain = self.domain.lower()
-        # Constant thicknesses
-        L_inner = self.phase_param.L_inner_0
-        L_outer = self.phase_param.L_outer_0
-        variables = self._get_standard_thickness_variables(L_inner, L_outer)
+        if self.reaction_loc == "interface":
+            c_inner = self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+            c_outer = self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+        else:
+            c_inner = self.phase_param.a_typ * (
+                self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+            )
+            c_outer = self.phase_param.a_typ * (
+                self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+            )
+        variables = self._get_standard_concentration_variables(c_inner, c_outer)
 
         # Reactions
         if self.reaction_loc == "interface":
@@ -47,12 +55,3 @@ def get_fundamental_variables(self):
         variables.update(self._get_standard_reaction_variables(zero, zero))
 
         return variables
-
-    def get_coupled_variables(self, variables):
-        # Concentrations (derived from thicknesses)
-        variables.update(self._get_standard_concentration_variables(variables))
-
-        # Add other standard coupled variables
-        variables.update(super().get_coupled_variables(variables))
-
-        return variables
diff --git a/pybamm/models/submodels/interface/sei/no_sei.py b/pybamm/models/submodels/interface/sei/no_sei.py
index 49d11a35e5..74c63c3abd 100644
--- a/pybamm/models/submodels/interface/sei/no_sei.py
+++ b/pybamm/models/submodels/interface/sei/no_sei.py
@@ -36,12 +36,6 @@ def get_fundamental_variables(self):
             zero = pybamm.FullBroadcast(
                 pybamm.Scalar(0), f"{domain} electrode", "current collector"
             )
-        variables = self._get_standard_thickness_variables(zero, zero)
+        variables = self._get_standard_concentration_variables(zero, zero)
         variables.update(self._get_standard_reaction_variables(zero, zero))
         return variables
-
-    def get_coupled_variables(self, variables):
-        variables.update(self._get_standard_concentration_variables(variables))
-        # Update whole cell variables, which also updates the "sum of" variables
-        variables.update(super().get_coupled_variables(variables))
-        return variables
diff --git a/pybamm/models/submodels/interface/sei/sei_growth.py b/pybamm/models/submodels/interface/sei/sei_growth.py
index 89bb662fbc..a25299f7d5 100644
--- a/pybamm/models/submodels/interface/sei/sei_growth.py
+++ b/pybamm/models/submodels/interface/sei/sei_growth.py
@@ -29,11 +29,14 @@ class SEIGrowth(BaseModel):
         Whether this is a submodel for standard SEI or SEI on cracks
     """
 
-    def __init__(
-        self, param, domain, reaction_loc, options, phase="primary", cracks=False
-    ):
+    def __init__(self, param, domain, options, phase="primary", cracks=False):
         super().__init__(param, domain, options=options, phase=phase, cracks=cracks)
-        self.reaction_loc = reaction_loc
+        if self.options.electrode_types[domain] == "planar":
+            self.reaction_loc = "interface"
+        elif self.options["x-average side reactions"] == "true":
+            self.reaction_loc = "x-average"
+        else:
+            self.reaction_loc = "full electrode"
         SEI_option = getattr(self.options, domain)["SEI"]
         if SEI_option == "ec reaction limited":
             pybamm.citations.register("Yang2017")
@@ -42,40 +45,52 @@ def __init__(
 
     def get_fundamental_variables(self):
         domain, Domain = self.domain_Domain
-        Ls = []
+        cs = []
         for pos in ["inner", "outer"]:
-            scale = self.phase_param.L_sei_0
             if self.reaction_loc == "x-average":
-                L_av = pybamm.Variable(
-                    f"X-averaged {domain} {pos} {self.reaction_name}thickness [m]",
+                c_sei_0 = self.phase_param.a_typ * (
+                    self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                    + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+                )
+                c_av = pybamm.Variable(
+                    f"X-averaged {domain} {pos} {self.reaction_name}"
+                    "concentration [mol.m-3]",
                     domain="current collector",
-                    scale=scale,
+                    scale=c_sei_0,
                 )
-                L_av.print_name = f"L_{pos}_av"
-                L = pybamm.PrimaryBroadcast(L_av, f"{domain} electrode")
+                c_av.print_name = f"c_{pos}_av"
+                c = pybamm.PrimaryBroadcast(c_av, f"{domain} electrode")
             elif self.reaction_loc == "full electrode":
-                L = pybamm.Variable(
-                    f"{Domain} {pos} {self.reaction_name}thickness [m]",
+                c_sei_0 = self.phase_param.a_typ * (
+                    self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                    + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+                )
+                c = pybamm.Variable(
+                    f"{Domain} {pos} {self.reaction_name}concentration [mol.m-3]",
                     domain=f"{domain} electrode",
                     auxiliary_domains={"secondary": "current collector"},
-                    scale=scale,
+                    scale=c_sei_0,
                 )
             elif self.reaction_loc == "interface":
-                L = pybamm.Variable(
-                    f"{Domain} {pos} {self.reaction_name}thickness [m]",
+                c_sei_0 = (
+                    self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                    + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+                )
+                c = pybamm.Variable(
+                    f"{Domain} {pos} {self.reaction_name}concentration [mol.m-2]",
                     domain="current collector",
-                    scale=scale,
+                    scale=c_sei_0,
                 )
-            L.print_name = f"L_{pos}"
-            Ls.append(L)
+            c.print_name = f"c_{pos}"
+            cs.append(c)
 
-        L_inner, L_outer = Ls
+        c_inner, c_outer = cs
 
         SEI_option = getattr(self.options, domain)["SEI"]
         if SEI_option.startswith("ec reaction limited"):
-            L_inner = 0 * L_inner  # Set L_inner to zero, copying domains
+            c_inner = 0 * c_inner  # Set L_inner to zero, copying domains
 
-        variables = self._get_standard_thickness_variables(L_inner, L_outer)
+        variables = self._get_standard_concentration_variables(c_inner, c_outer)
 
         return variables
 
@@ -182,7 +197,6 @@ def get_coupled_variables(self, variables):
         j_inner = inner_sei_proportion * Arrhenius * j_sei
         j_outer = (1 - inner_sei_proportion) * Arrhenius * j_sei
 
-        variables.update(self._get_standard_concentration_variables(variables))
         variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
 
         # Add other standard coupled variables
@@ -195,12 +209,28 @@ def set_rhs(self, variables):
         param = self.param
         domain, Domain = self.domain_Domain
 
-        if self.reaction_loc == "x-average":
-            L_inner = variables[
-                f"X-averaged {domain} inner {self.reaction_name}thickness [m]"
+        if self.reaction_loc == "interface":
+            c_inner = variables[
+                f"{Domain} inner {self.reaction_name}concentration [mol.m-2]"
+            ]
+            c_outer = variables[
+                f"{Domain} outer {self.reaction_name}concentration [mol.m-2]"
+            ]
+            j_inner = variables[
+                f"{Domain} electrode inner {self.reaction_name}"
+                "interfacial current density [A.m-2]"
+            ]
+            j_outer = variables[
+                f"{Domain} electrode outer {self.reaction_name}"
+                "interfacial current density [A.m-2]"
             ]
-            L_outer = variables[
-                f"X-averaged {domain} outer {self.reaction_name}thickness [m]"
+            a = 1
+        elif self.reaction_loc == "x-average":
+            c_inner = variables[
+                f"X-averaged {domain} inner {self.reaction_name}concentration [mol.m-3]"
+            ]
+            c_outer = variables[
+                f"X-averaged {domain} outer {self.reaction_name}concentration [mol.m-3]"
             ]
             j_inner = variables[
                 f"X-averaged {domain} electrode inner {self.reaction_name}"
@@ -210,10 +240,17 @@ def set_rhs(self, variables):
                 f"X-averaged {domain} electrode outer {self.reaction_name}"
                 "interfacial current density [A.m-2]"
             ]
-
+            a = variables[
+                f"X-averaged {domain} electrode {self.phase_name}"
+                "surface area to volume ratio [m-1]"
+            ]
         else:
-            L_inner = variables[f"{Domain} inner {self.reaction_name}thickness [m]"]
-            L_outer = variables[f"{Domain} outer {self.reaction_name}thickness [m]"]
+            c_inner = variables[
+                f"{Domain} inner {self.reaction_name}concentration [mol.m-3]"
+            ]
+            c_outer = variables[
+                f"{Domain} outer {self.reaction_name}concentration [mol.m-3]"
+            ]
             j_inner = variables[
                 f"{Domain} electrode inner {self.reaction_name}"
                 "interfacial current density [A.m-2]"
@@ -222,70 +259,73 @@ def set_rhs(self, variables):
                 f"{Domain} electrode outer {self.reaction_name}"
                 "interfacial current density [A.m-2]"
             ]
+            a = variables[
+                f"{Domain} electrode {self.phase_name}"
+                "surface area to volume ratio [m-1]"
+            ]
 
-        # The spreading term acts to spread out SEI along the cracks as they grow.
-        # For SEI on initial surface (as opposed to cracks), it is zero.
         if self.reaction == "SEI on cracks":
             if self.reaction_loc == "x-average":
-                l_cr = variables[f"X-averaged {domain} particle crack length [m]"]
-                dl_cr = variables[f"X-averaged {domain} particle cracking rate [m.s-1]"]
+                roughness = variables[f"X-averaged {domain} electrode roughness ratio"]
             else:
-                l_cr = variables[f"{Domain} particle crack length [m]"]
-                dl_cr = variables[f"{Domain} particle cracking rate [m.s-1]"]
-            spreading_outer = (
-                dl_cr / l_cr * (self.phase_param.L_outer_crack_0 - L_outer)
-            )
-            spreading_inner = (
-                dl_cr / l_cr * (self.phase_param.L_inner_crack_0 - L_inner)
-            )
-        else:
-            spreading_outer = 0
-            spreading_inner = 0
+                roughness = variables[f"{Domain} electrode roughness ratio"]
+            a *= roughness - 1  # Replace surface area with crack area
 
         # a * j_sei / F is the rate of consumption of li moles by SEI reaction
         # 1/z_sei converts from li moles to SEI moles (z_sei=li mol per sei mol)
         # a * j_sei / (F * z_sei) is the rate of consumption of SEI moles by SEI
         # reaction
-        # V_bar / a converts from SEI moles to SEI thickness
-        # V_bar * j_sei / (F * z_sei) is the rate of SEI thickness change
-        dLdt_SEI_inner = (
-            phase_param.V_bar_inner * j_inner / (param.F * phase_param.z_sei)
-        )
-        dLdt_SEI_outer = (
-            phase_param.V_bar_outer * j_outer / (param.F * phase_param.z_sei)
-        )
-
-        # we have to add the spreading rate to account for cracking
+        dcdt_SEI_inner = a * j_inner / (param.F * phase_param.z_sei)
+        dcdt_SEI_outer = a * j_outer / (param.F * phase_param.z_sei)
+
         SEI_option = getattr(self.options, domain)["SEI"]
         if SEI_option.startswith("ec reaction limited"):
-            self.rhs = {L_outer: -dLdt_SEI_outer + spreading_outer}
+            self.rhs = {c_outer: -dcdt_SEI_outer}
         else:
-            self.rhs = {
-                L_inner: -dLdt_SEI_inner + spreading_inner,
-                L_outer: -dLdt_SEI_outer + spreading_outer,
-            }
+            self.rhs = {c_inner: -dcdt_SEI_inner, c_outer: -dcdt_SEI_outer}
 
     def set_initial_conditions(self, variables):
         domain, Domain = self.domain_Domain
-        if self.reaction_loc == "x-average":
-            L_inner = variables[
-                f"X-averaged {domain} inner {self.reaction_name}thickness [m]"
+        if self.reaction_loc == "interface":
+            c_inner = variables[
+                f"{Domain} inner {self.reaction_name}concentration [mol.m-2]"
+            ]
+            c_outer = variables[
+                f"{Domain} outer {self.reaction_name}concentration [mol.m-2]"
             ]
-            L_outer = variables[
-                f"X-averaged {domain} outer {self.reaction_name}thickness [m]"
+            a = 1
+        elif self.reaction_loc == "x-average":
+            c_inner = variables[
+                f"X-averaged {domain} inner {self.reaction_name}concentration [mol.m-3]"
             ]
+            c_outer = variables[
+                f"X-averaged {domain} outer {self.reaction_name}concentration [mol.m-3]"
+            ]
+            a = self.phase_param.a_typ
         else:
-            L_inner = variables[f"{Domain} inner {self.reaction_name}thickness [m]"]
-            L_outer = variables[f"{Domain} outer {self.reaction_name}thickness [m]"]
+            c_inner = variables[
+                f"{Domain} inner {self.reaction_name}concentration [mol.m-3]"
+            ]
+            c_outer = variables[
+                f"{Domain} outer {self.reaction_name}concentration [mol.m-3]"
+            ]
+            a = self.phase_param.a_typ
 
         if self.reaction == "SEI on cracks":
             L_inner_0 = self.phase_param.L_inner_crack_0
             L_outer_0 = self.phase_param.L_outer_crack_0
+            roughness_init = 1 + 2 * (
+                self.param.n.l_cr_0 * self.param.n.rho_cr * self.param.n.w_cr
+            )
+            a *= roughness_init - 1  # Replace surface area with crack area
         else:
             L_inner_0 = self.phase_param.L_inner_0
             L_outer_0 = self.phase_param.L_outer_0
+
+        c_inner_0 = a * L_inner_0 / self.phase_param.V_bar_inner
+        c_outer_0 = a * L_outer_0 / self.phase_param.V_bar_outer
         SEI_option = getattr(self.options, domain)["SEI"]
         if SEI_option.startswith("ec reaction limited"):
-            self.initial_conditions = {L_outer: L_inner_0 + L_outer_0}
+            self.initial_conditions = {c_outer: c_inner_0 + c_outer_0}
         else:
-            self.initial_conditions = {L_inner: L_inner_0, L_outer: L_outer_0}
+            self.initial_conditions = {c_inner: c_inner_0, c_outer: c_outer_0}
diff --git a/pybamm/models/submodels/interface/sei/sei_thickness.py b/pybamm/models/submodels/interface/sei/sei_thickness.py
new file mode 100644
index 0000000000..699ee253f6
--- /dev/null
+++ b/pybamm/models/submodels/interface/sei/sei_thickness.py
@@ -0,0 +1,200 @@
+#
+# Class for converting SEI concentration into thickness
+#
+import pybamm
+from .base_sei import BaseModel
+
+
+class SEIThickness(BaseModel):
+    """
+    Class for converting SEI concentration into thickness
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    reaction_loc : str
+        Where the reaction happens: "x-average" (SPM, SPMe, etc),
+        "full electrode" (full DFN), or "interface" (half-cell model)
+    options : dict
+        A dictionary of options to be passed to the model.
+    phase : str, optional
+        Phase of the particle (default is "primary")
+    cracks : bool, optional
+        Whether this is a submodel for standard SEI or SEI on cracks
+    """
+
+    def __init__(self, param, domain, options, phase="primary", cracks=False):
+        super().__init__(param, domain, options=options, phase=phase, cracks=cracks)
+        if self.options.electrode_types[domain] == "planar":
+            self.reaction_loc = "interface"
+        elif self.options["x-average side reactions"] == "true":
+            self.reaction_loc = "x-average"
+        else:
+            self.reaction_loc = "full electrode"
+
+    def get_coupled_variables(self, variables):
+        """Update variables related to the SEI thickness."""
+        domain, Domain = self.domain_Domain
+        phase_param = self.phase_param
+        reaction_name = self.reaction_name
+        SEI_option = getattr(self.options, domain)["SEI"]
+        crack_option = getattr(self.options, domain)["SEI on cracks"]
+        if self.options["working electrode"] != "both" and domain == "negative":
+            crack_option = "false"  # required if SEI on cracks is used for half-cells
+
+        # Set scales to one for the "no SEI" model so that they are not required
+        # by parameter values in general
+        if SEI_option == "none":
+            c_to_L_inner = 1
+            c_to_L_outer = 1
+            R_sei = 1
+        else:
+            if self.reaction_loc == "interface":
+                c_to_L_inner = phase_param.V_bar_inner
+                c_to_L_outer = phase_param.V_bar_outer
+            else:
+                a = variables[
+                    f"{Domain} electrode {self.phase_name}"
+                    "surface area to volume ratio [m-1]"
+                ]
+                c_to_L_inner = phase_param.V_bar_inner / a
+                c_to_L_outer = phase_param.V_bar_outer / a
+            R_sei = phase_param.R_sei
+
+        if self.reaction_loc == "interface":
+            # c is an interfacial quantity [mol.m-2]
+            c_inner = variables[
+                f"{Domain} inner {reaction_name}concentration [mol.m-2]"
+            ]
+            c_outer = variables[
+                f"{Domain} outer {reaction_name}concentration [mol.m-2]"
+            ]
+        else:
+            # c is a bulk quantity [mol.m-3]
+            c_inner = variables[
+                f"{Domain} inner {reaction_name}concentration [mol.m-3]"
+            ]
+            c_outer = variables[
+                f"{Domain} outer {reaction_name}concentration [mol.m-3]"
+            ]
+
+        if self.reaction == "SEI":
+            L_inner = c_inner * c_to_L_inner  # inner SEI thickness
+            L_outer = c_outer * c_to_L_outer  # outer SEI thickness
+
+            L_inner_av = pybamm.x_average(L_inner)
+            L_outer_av = pybamm.x_average(L_outer)
+
+            L_SEI = L_inner + L_outer  # SEI thickness
+            L_SEI_av = pybamm.x_average(L_SEI)
+
+            variables.update(
+                {
+                    f"X-averaged {self.domain} electrode resistance "
+                    "[Ohm.m2]": L_SEI_av * R_sei,
+                }
+            )
+        # Thickness variables are handled slightly differently for SEI on cracks
+        elif self.reaction == "SEI on cracks":
+            # if SEI on cracks is false, skip over roughness to avoid division by zero
+            if crack_option == "false":
+                L_inner = c_inner * c_to_L_inner
+                L_outer = c_outer * c_to_L_outer
+            else:
+                roughness = variables[f"{Domain} electrode roughness ratio"]
+                L_inner = c_inner * c_to_L_inner / (roughness - 1)
+                L_outer = c_outer * c_to_L_outer / (roughness - 1)
+
+            L_inner_av = pybamm.x_average(L_inner)
+            L_outer_av = pybamm.x_average(L_outer)
+
+            L_SEI = L_inner + L_outer
+            L_SEI_av = pybamm.x_average(L_SEI)
+
+        variables.update(
+            {
+                f"{Domain} inner {reaction_name}thickness [m]": L_inner,
+                f"X-averaged {domain} inner {reaction_name}thickness [m]": L_inner_av,
+                f"{Domain} outer {reaction_name}thickness [m]": L_outer,
+                f"X-averaged {domain} outer {reaction_name}thickness [m]": L_outer_av,
+                f"{Domain} {reaction_name}thickness [m]": L_SEI,
+                f"X-averaged {domain} {reaction_name}thickness [m]": L_SEI_av,
+                f"{Domain} total {reaction_name}thickness [m]": L_SEI,
+                f"X-averaged {domain} total {reaction_name}thickness [m]": L_SEI_av,
+            }
+        )
+
+        # Calculate change in total SEI moles with respect to initial state
+        # If there is no SEI, skip and return 0 because parameters may be undefined
+        if self.reaction == "SEI" and SEI_option == "none":
+            Q_sei = pybamm.Scalar(0)
+        elif self.reaction == "SEI on cracks" and crack_option == "false":
+            Q_sei = pybamm.Scalar(0)
+        else:
+            if self.reaction_loc == "interface":
+                # c is an interfacial quantity [mol.m-2]
+                c_sei = variables[
+                    f"{Domain} total {self.reaction_name}concentration [mol.m-2]"
+                ]
+                c_sei_av = pybamm.yz_average(c_sei)
+                c_sei_0 = (
+                    self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                    + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+                )
+            else:
+                # c is a bulk quantity [mol.m-3]
+                c_sei = variables[
+                    f"X-averaged {domain} total {self.reaction_name}"
+                    "concentration [mol.m-3]"
+                ]
+                c_sei_av = pybamm.yz_average(c_sei)
+                c_sei_0 = self.phase_param.a_typ * (
+                    self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                    + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+                )
+            z_sei = self.phase_param.z_sei
+            if self.reaction == "SEI":
+                delta_c_SEI = c_sei_av - c_sei_0
+            elif self.reaction == "SEI on cracks":
+                if domain == "negative":
+                    roughness_init = 1 + 2 * (
+                        self.param.n.l_cr_0 * self.param.n.rho_cr * self.param.n.w_cr
+                    )
+                elif domain == "positive":
+                    roughness_init = 1 + 2 * (
+                        self.param.p.l_cr_0 * self.param.p.rho_cr * self.param.p.w_cr
+                    )
+                c_sei_cracks_0 = (
+                    (roughness_init - 1)
+                    * self.phase_param.a_typ
+                    * (
+                        self.phase_param.L_inner_crack_0 / self.phase_param.V_bar_inner
+                        + self.phase_param.L_outer_crack_0
+                        / self.phase_param.V_bar_outer
+                    )
+                )
+                delta_c_SEI = c_sei_av - c_sei_cracks_0
+
+            if self.reaction_loc == "interface":
+                L_k = 1
+            elif domain == "negative":
+                L_k = self.param.n.L
+            elif domain == "positive":
+                L_k = self.param.p.L
+
+            # Multiply delta_n_SEI by V_k to get total moles of SEI formed
+            # Multiply by z_sei to get total lithium moles consumed by SEI
+            V_k = L_k * self.param.L_y * self.param.L_z
+            Q_sei = delta_c_SEI * V_k * z_sei
+
+        variables.update(
+            {
+                f"Loss of lithium to {domain} {self.reaction_name}[mol]": Q_sei,
+                f"Loss of capacity to {domain} {self.reaction_name}[A.h]": Q_sei
+                * self.param.F
+                / 3600,
+            }
+        )
+
+        return variables
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 083e3b648b..7eed921d2d 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
@@ -34,6 +34,43 @@ def positive_radius(x):
         param["Current function [A]"] = 0.5 * param["Nominal cell capacity [A.h]"]
         self.run_basic_processing_test({}, parameter_values=param)
 
+    def test_conservation_coupled_degradation(self):
+        # Test that lithium is conserved when all degradation mechanisms are enabled
+        model = pybamm.lithium_ion.DFN(
+            {
+                "SEI": "reaction limited",
+                "SEI porosity change": "true",
+                "lithium plating": "partially reversible",
+                "particle mechanics": ("swelling and cracking", "swelling only"),
+                "SEI on cracks": "true",
+                "loss of active material": "stress and reaction-driven",
+            }
+        )
+        param = pybamm.ParameterValues("OKane2022")
+
+        # optimize mesh size for mechanical degradation
+        var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_n": 30, "r_p": 20}
+        solver = pybamm.CasadiSolver(mode="fast")
+
+        # solve
+        sim = pybamm.Simulation(
+            model, parameter_values=param, var_pts=var_pts, solver=solver
+        )
+        solution = sim.solve([0, 3500])
+        BOD = solution["Total lithium in particles [mol]"].entries[0]
+        particles = solution["Total lithium in particles [mol]"].entries[-1]
+        neg = solution[
+            "Loss of lithium due to loss of active material in negative electrode [mol]"
+        ].entries[-1]
+        pos = solution[
+            "Loss of lithium due to loss of active material in positive electrode [mol]"
+        ].entries[-1]
+        side = solution["Total lithium lost to side reactions [mol]"].entries[-1]
+        EOD = particles + neg + pos + side
+
+        # compare
+        np.testing.assert_array_almost_equal(BOD, EOD, decimal=8)
+
 
 class TestDFNWithSizeDistribution(TestCase):
     def setUp(self):