diff --git a/docs/source/examples/notebooks/models/thermal-models.ipynb b/docs/source/examples/notebooks/models/thermal-models.ipynb
index 788a02a88f..9ee4c1328a 100644
--- a/docs/source/examples/notebooks/models/thermal-models.ipynb
+++ b/docs/source/examples/notebooks/models/thermal-models.ipynb
@@ -519,7 +519,7 @@
   ],
   "metadata": {
     "kernelspec": {
-      "display_name": "Python 3 (ipykernel)",
+      "display_name": "dev",
       "language": "python",
       "name": "python3"
     },
@@ -533,7 +533,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.9.0"
+      "version": "3.9.16"
     },
     "toc": {
       "base_numbering": 1,
@@ -547,6 +547,11 @@
       "toc_position": {},
       "toc_section_display": true,
       "toc_window_display": true
+    },
+    "vscode": {
+      "interpreter": {
+        "hash": "bca2b99bfac80e18288b793d52fa0653ab9b5fe5d22e7b211c44eb982a41c00c"
+      }
     }
   },
   "nbformat": 4,
diff --git a/examples/scripts/thermal_lithium_ion.py b/examples/scripts/thermal_lithium_ion.py
index 9240c6e895..cd950cd701 100644
--- a/examples/scripts/thermal_lithium_ion.py
+++ b/examples/scripts/thermal_lithium_ion.py
@@ -3,58 +3,35 @@
 #
 
 import pybamm
-import numpy as np
 
-# load model
 pybamm.set_logging_level("INFO")
 
-options = {"thermal": "x-full"}
-full_thermal_model = pybamm.lithium_ion.SPMe(options)
-
-options = {"thermal": "x-lumped"}
-lumped_thermal_model = pybamm.lithium_ion.SPMe(options)
-
-models = [full_thermal_model, lumped_thermal_model]
-
-# load parameter values and process models and geometry
-param = models[0].default_parameter_values
+# load models
+models = [
+    pybamm.lithium_ion.SPMe({"thermal": "x-full"}),
+    pybamm.lithium_ion.SPMe({"thermal": "x-lumped"}),
+]
 
-# for x-full, cooling is only implemented on the surfaces
-# so set other forms of cooling to zero for comparison.
-param.update(
+# load parameter values and update cooling coefficients
+parameter_values = pybamm.ParameterValues("Marquis2019")
+parameter_values.update(
     {
-        "Negative current collector"
-        + " surface heat transfer coefficient [W.m-2.K-1]": 5,
-        "Positive current collector"
-        + " surface heat transfer coefficient [W.m-2.K-1]": 5,
-        "Negative tab heat transfer coefficient [W.m-2.K-1]": 0,
-        "Positive tab heat transfer coefficient [W.m-2.K-1]": 0,
-        "Edge heat transfer coefficient [W.m-2.K-1]": 0,
+        "Negative current collector surface heat transfer coefficient [W.m-2.K-1]"
+        "": 5,
+        "Positive current collector surface heat transfer coefficient [W.m-2.K-1]"
+        "": 5,
+        "Negative tab heat transfer coefficient [W.m-2.K-1]": 10,
+        "Positive tab heat transfer coefficient [W.m-2.K-1]": 10,
+        "Edge heat transfer coefficient [W.m-2.K-1]": 5,
     }
 )
 
+# create and solve simulations
+sols = []
 for model in models:
-    param.process_model(model)
-
-# set mesh
-var_pts = {"x_n": 10, "x_s": 10, "x_p": 10, "r_n": 10, "r_p": 10}
-
-# discretise models
-for model in models:
-    # create geometry
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve model
-solutions = [None] * len(models)
-t_eval = np.linspace(0, 3600, 100)
-for i, model in enumerate(models):
-    solver = pybamm.ScipySolver(atol=1e-8, rtol=1e-8)
-    solution = solver.solve(model, t_eval)
-    solutions[i] = solution
+    sim = pybamm.Simulation(model, parameter_values=parameter_values)
+    sol = sim.solve([0, 3600])
+    sols.append(sol)
 
 # plot
 output_variables = [
@@ -63,5 +40,4 @@
     "Cell temperature [K]",
 ]
 labels = ["Full thermal model", "Lumped thermal model"]
-plot = pybamm.QuickPlot(solutions, output_variables, labels)
-plot.dynamic_plot()
+pybamm.dynamic_plot(sols, output_variables, labels=labels)
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index db0e38710a..bf0600577a 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -26,8 +26,8 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "true" or "false". "false" is the default, since calculating discharge
                 energy can be computationally expensive for simple models like SPM.
             * "cell geometry" : str
-                Sets the geometry of the cell. Can be "pouch" (default) or
-                "arbitrary". The arbitrary geometry option solves a 1D electrochemical
+                Sets the geometry of the cell. Can be "arbitrary" (default) or
+                "pouch". The arbitrary geometry option solves a 1D electrochemical
                 model with prescribed cell volume and cross-sectional area, and
                 (if thermal effects are included) solves a lumped thermal model
                 with prescribed surface area for cooling.
@@ -283,9 +283,9 @@ def __init__(self, extra_options):
         default_options = {
             name: options[0] for name, options in self.possible_options.items()
         }
+        extra_options = extra_options or {}
 
         # Change the default for cell geometry based on which thermal option is provided
-        extra_options = extra_options or {}
         # return "none" if option not given
         thermal_option = extra_options.get("thermal", "none")
         if thermal_option in ["none", "isothermal", "lumped"]:
diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
index 54c4a1b467..2ac83f1aef 100644
--- a/pybamm/models/submodels/thermal/base_thermal.py
+++ b/pybamm/models/submodels/thermal/base_thermal.py
@@ -175,6 +175,8 @@ def _get_standard_coupled_variables(self, variables):
                 "Total heating [W.m-3]": Q,
                 "X-averaged total heating [W.m-3]": Q_av,
                 "Volume-averaged total heating [W.m-3]": Q_vol_av,
+                "Negative current collector Ohmic heating [W.m-3]": Q_ohm_s_cn,
+                "Positive current collector Ohmic heating [W.m-3]": Q_ohm_s_cp,
             }
         )
         return variables
diff --git a/pybamm/models/submodels/thermal/lumped.py b/pybamm/models/submodels/thermal/lumped.py
index 0b7387390f..18d5d1bea2 100644
--- a/pybamm/models/submodels/thermal/lumped.py
+++ b/pybamm/models/submodels/thermal/lumped.py
@@ -54,7 +54,6 @@ def set_rhs(self, variables):
         T_amb = variables["Ambient temperature [K]"]
 
         # Account for surface area to volume ratio in cooling coefficient
-        # The factor 1/delta^2 comes from the choice of non-dimensionalisation.
         if self.options["cell geometry"] == "pouch":
             cell_volume = self.param.L * self.param.L_y * self.param.L_z
 
diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
index 35d226d6e3..4117ae84f1 100644
--- a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
+++ b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
@@ -57,8 +57,7 @@ def set_rhs(self, variables):
         T_amb = variables["Ambient temperature [K]"]
 
         # Account for surface area to volume ratio of pouch cell in cooling
-        # coefficient. Note: the factor 1/delta^2 comes from the choice of
-        # non-dimensionalisation
+        # coefficient.
         yz_surface_area = self.param.L_y * self.param.L_z
         cell_volume = self.param.L * self.param.L_y * self.param.L_z
         yz_surface_cooling_coefficient = (
diff --git a/pybamm/models/submodels/thermal/x_full.py b/pybamm/models/submodels/thermal/x_full.py
index 49d83a6f2e..da95e2a46a 100644
--- a/pybamm/models/submodels/thermal/x_full.py
+++ b/pybamm/models/submodels/thermal/x_full.py
@@ -41,8 +41,16 @@ def get_fundamental_variables(self):
             T_dict[domain] = T_k
 
         T = pybamm.concatenation(*T_dict.values())
-        T_cn = pybamm.boundary_value(T_dict["negative electrode"], "left")
-        T_cp = pybamm.boundary_value(T_dict["positive electrode"], "right")
+        T_cn = pybamm.Variable(
+            "Negative current collector temperature [K]",
+            domain="current collector",
+            scale=self.param.T_ref,
+        )
+        T_cp = pybamm.Variable(
+            "Positive current collector temperature [K]",
+            domain="current collector",
+            scale=self.param.T_ref,
+        )
         T_x_av = self._x_average(T, T_cn, T_cp)
         T_vol_av = self._yz_average(T_x_av)
         T_dict.update(
@@ -63,57 +71,110 @@ def get_coupled_variables(self, variables):
 
     def set_rhs(self, variables):
         T = variables["Cell temperature [K]"]
+        T_cn = variables["Negative current collector temperature [K]"]
         T_n = variables["Negative electrode temperature [K]"]
         T_s = variables["Separator temperature [K]"]
         T_p = variables["Positive electrode temperature [K]"]
-
+        T_cp = variables["Positive current collector temperature [K]"]
         Q = variables["Total heating [W.m-3]"]
+        Q_cn = variables["Negative current collector Ohmic heating [W.m-3]"]
+        Q_cp = variables["Positive current collector Ohmic heating [W.m-3]"]
+        T_amb = variables["Ambient temperature [K]"]
 
-        # Define volumetric heat capacity
+        # Define volumetric heat capacity for electrode/separator/electrode sandwich
         rho_c_p = pybamm.concatenation(
             self.param.n.rho_c_p(T_n),
             self.param.s.rho_c_p(T_s),
             self.param.p.rho_c_p(T_p),
         )
 
-        # Define thermal conductivity
+        # Define thermal conductivity for electrode/separator/electrode sandwich
         lambda_ = pybamm.concatenation(
             self.param.n.lambda_(T_n),
             self.param.s.lambda_(T_s),
             self.param.p.lambda_(T_p),
         )
 
-        # Fourier's law for heat flux
-        q = -lambda_ * pybamm.grad(T)
+        # Calculate edge/tab cooling
+        L_y = self.param.L_y
+        L_z = self.param.L_z
+        L_cn = self.param.n.L_cc
+        L_cp = self.param.p.L_cc
+        h_cn = self.param.n.h_cc
+        h_cp = self.param.p.h_cc
+        lambda_n = self.param.n.lambda_(T_n)
+        lambda_p = self.param.p.lambda_(T_p)
+        # Negative current collector
+        volume_cn = L_cn * L_y * L_z
+        negative_tab_area = self.param.n.L_tab * self.param.n.L_cc
+        edge_area_cn = 2 * (L_y + L_z) * L_cn - negative_tab_area
+        negative_tab_cooling_coefficient = (
+            -self.param.n.h_tab * negative_tab_area / volume_cn
+        )
+        edge_cooling_coefficient_cn = -self.param.h_edge * edge_area_cn / volume_cn
+        cooling_coefficient_cn = (
+            negative_tab_cooling_coefficient + edge_cooling_coefficient_cn
+        )
+        # Electrode/separator/electrode sandwich
+        area_to_volume = (
+            2 * (self.param.L_y + self.param.L_z) / (self.param.L_y * self.param.L_z)
+        )
+        cooling_coefficient = -self.param.h_edge * area_to_volume
+        # Positive current collector
+        volume_cp = L_cp * L_y * L_z
+        positive_tab_area = self.param.p.L_tab * self.param.p.L_cc
+        edge_area_cp = 2 * (L_y + L_z) * L_cp - positive_tab_area
+        positive_tab_cooling_coefficient = (
+            -self.param.p.h_tab * positive_tab_area / volume_cp
+        )
+        edge_cooling_coefficient_cp = -self.param.h_edge * edge_area_cp / volume_cp
+        cooling_coefficient_cp = (
+            positive_tab_cooling_coefficient + edge_cooling_coefficient_cp
+        )
 
-        # N.B only y-z surface cooling is implemented for this model
-        self.rhs = {T: (-pybamm.div(q) + Q) / rho_c_p}
+        self.rhs = {
+            T_cn: (
+                (
+                    pybamm.boundary_value(lambda_n, "left")
+                    * pybamm.boundary_gradient(T_n, "left")
+                    - h_cn * (T_cn - T_amb)
+                )
+                / L_cn
+                + Q_cn
+                + cooling_coefficient_cn * (T_cn - T_amb)
+            )
+            / self.param.n.rho_c_p_cc(T_cn),
+            T: (
+                pybamm.div(lambda_ * pybamm.grad(T))
+                + Q
+                + cooling_coefficient * (T - T_amb)
+            )
+            / rho_c_p,
+            T_cp: (
+                (
+                    -pybamm.boundary_value(lambda_p, "right")
+                    * pybamm.boundary_gradient(T_p, "right")
+                    - h_cp * (T_cp - T_amb)
+                )
+                / L_cp
+                + Q_cp
+                + cooling_coefficient_cp * (T_cp - T_amb)
+            )
+            / self.param.p.rho_c_p_cc(T_cp),
+        }
 
     def set_boundary_conditions(self, variables):
         T = variables["Cell temperature [K]"]
-        T_n_left = pybamm.boundary_value(T, "left")
-        T_p_right = pybamm.boundary_value(T, "right")
-        T_amb = variables["Ambient temperature [K]"]
+        T_cn = variables["Negative current collector temperature [K]"]
+        T_cp = variables["Positive current collector temperature [K]"]
 
-        # N.B only y-z surface cooling is implemented for this thermal model.
-        # Tab and edge cooling is not accounted for.
         self.boundary_conditions = {
-            T: {
-                "left": (
-                    self.param.n.h_cc
-                    * (T_n_left - T_amb)
-                    / self.param.n.lambda_(T_n_left),
-                    "Neumann",
-                ),
-                "right": (
-                    -self.param.p.h_cc
-                    * (T_p_right - T_amb)
-                    / self.param.p.lambda_(T_p_right),
-                    "Neumann",
-                ),
-            }
+            T: {"left": (T_cn, "Dirichlet"), "right": (T_cp, "Dirichlet")}
         }
 
     def set_initial_conditions(self, variables):
         T = variables["Cell temperature [K]"]
-        self.initial_conditions = {T: self.param.T_init}
+        T_cn = variables["Negative current collector temperature [K]"]
+        T_cp = variables["Positive current collector temperature [K]"]
+        T_init = self.param.T_init
+        self.initial_conditions = {T_cn: T_init, T: T_init, T_cp: T_init}
diff --git a/pybamm/parameters/lead_acid_parameters.py b/pybamm/parameters/lead_acid_parameters.py
index b1302aff05..7da6a19410 100644
--- a/pybamm/parameters/lead_acid_parameters.py
+++ b/pybamm/parameters/lead_acid_parameters.py
@@ -249,9 +249,6 @@ def _set_parameters(self):
 
         # Macroscale geometry
         self.L = self.geo.L
-        # In lead-acid the current collector and electrodes are the same (same
-        # thickness)
-        self.L_cc = self.L
 
         # Thermal
         self.rho_c_p = self.therm.rho_c_p
@@ -265,10 +262,20 @@ def _set_parameters(self):
             self.b_e = self.geo.b_e
             self.epsilon_inactive = pybamm.Scalar(0)
             return
+
+        # In lead-acid the current collector and electrodes are the same (same
+        # thickness)
+        self.L_cc = self.L
         # for lead-acid the electrodes and current collector are the same
         self.rho_c_p_cc = self.therm.rho_c_p
         self.lambda_cc = self.therm.lambda_
 
+        # Tab geometry (for pouch cells)
+        self.L_tab = self.geo.L_tab
+        self.centre_y_tab = self.geo.centre_y_tab
+        self.centre_z_tab = self.geo.centre_z_tab
+        self.A_tab = self.geo.A_tab
+
         # Microstructure
         self.b_e = self.geo.b_e
         self.b_s = self.geo.b_s