diff --git a/examples/scripts/compare_lithium_ion_heat_of_mixing.py b/examples/scripts/compare_lithium_ion_heat_of_mixing.py
new file mode 100644
index 0000000000..0c65ad640d
--- /dev/null
+++ b/examples/scripts/compare_lithium_ion_heat_of_mixing.py
@@ -0,0 +1,43 @@
+#
+# Compare lithium-ion battery models
+#
+import pybamm
+
+pybamm.set_logging_level("INFO")
+# load models
+models = [
+    pybamm.lithium_ion.DFN({"heat of mixing": "true","thermal": "x-lumped"},name="hom"),
+    pybamm.lithium_ion.DFN({"thermal": "x-lumped"},name="nhom")
+]
+
+parameter_values = pybamm.ParameterValues("Chen2020_composite")
+
+# create and run simulations
+sims = []
+for model in models:
+    sim = pybamm.Simulation(model, parameter_values=parameter_values)
+    sim.solve([0, 4500])
+    sims.append(sim)
+
+# plot
+pybamm.dynamic_plot(
+    sims,
+    [
+        [
+            "Average negative primary particle concentration",
+            "Average negative secondary particle concentration",
+        ],
+        [
+            "X-averaged negative electrode primary volumetric "
+            "interfacial current density [A.m-3]",
+            "X-averaged negative electrode secondary volumetric "
+            "interfacial current density [A.m-3]",
+            "X-averaged negative electrode volumetric "
+            "interfacial current density [A.m-3]",
+        ],
+        "X-averaged negative electrode primary open-circuit potential [V]",
+        "X-averaged negative electrode secondary open-circuit potential [V]",
+        "Average positive particle concentration [mol.m-3]",
+        "Voltage [V]",
+    ],
+)
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index e82f3d5868..03246b4478 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -199,6 +199,7 @@ def __init__(self, extra_options):
                 "composite",
                 "integrated",
             ],
+            "heat of mixing": ["false", "true"],
             "hydrolysis": ["false", "true"],
             "intercalation kinetics": [
                 "symmetric Butler-Volmer",
@@ -414,6 +415,11 @@ def __init__(self, extra_options):
 
         # Options not yet compatible with particle-size distributions
         if options["particle size"] == "distribution":
+            if options["heat of mixing"] != "false":
+                raise NotImplementedError(
+                    "Heat of mixing submodels do not yet support particle-size "
+                    "distributions."
+                )
             if options["lithium plating"] != "none":
                 raise NotImplementedError(
                     "Lithium plating submodels do not yet support particle-size "
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index 2068167e8f..01ff43afe2 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -22,6 +22,7 @@
 'current collector': 'uniform' (possible: ['uniform', 'potential pair', 'potential pair quite conductive'])
 'dimensionality': 0 (possible: [0, 1, 2])
 'electrolyte conductivity': 'default' (possible: ['default', 'full', 'leading order', 'composite', 'integrated'])
+'heat of mixing': 'false' (possible: ['false', 'true'])
 'hydrolysis': 'false' (possible: ['false', 'true'])
 'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey'])
 'interface utilisation': 'full' (possible: ['full', 'constant', 'current-driven'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 99b2b10148..5e2c084e3a 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -67,6 +67,76 @@ def test_well_posed_thermal_2plus1D(self):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_isothermal_heat_source_hom(self):
+        options = {
+            "calculate heat source for isothermal models": "true",
+            "thermal": "isothermal",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_2plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 1,
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_model_1D_hom(self):
+        options = {"thermal": "x-lumped", "heat of mixing": "true"}
+        self.check_well_posedness(options)
+
+    def test_well_posed_x_full_thermal_model_hom(self):
+        options = {
+            "thermal": "x-full",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_1plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 1,
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_2plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_thermal_1plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 1,
+            "thermal": "x-lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_thermal_2plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "thermal": "x-lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_contact_resistance(self):
         options = {"contact resistance": "true"}
         self.check_well_posedness(options)