diff --git a/CHANGELOG.md b/CHANGELOG.md
index f5f90b01a2..52450aac82 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,11 @@
 ## Features
 
 -   Added "Discharge energy [W.h]", which is the integral of the power in Watts, as an optional output. Set the option "calculate discharge energy" to "true" to get this output ("false" by default, since it can slow down some of the simple models) ([#1969](https://github.com/pybamm-team/PyBaMM/pull/1969)))
+-   Added an option "calculate heat source for isothermal models" to choose whether or not the heat generation terms are computed when running models with the option `thermal="isothermal"`  ([#1958](https://github.com/pybamm-team/PyBaMM/pull/1958))
+
+## Bug fixes
+
+-   Fixed a bug where isothermal models did not compute any heat source terms ([#1958](https://github.com/pybamm-team/PyBaMM/pull/1958))
 
 ## Breaking changes
 
@@ -12,7 +17,7 @@
 
 ## Features
 
--   Isothermal models now compute heat source terms (but the temperature remains constant). The models now also account for current collector heating when `dimensionality=0` ([#1929](https://github.com/pybamm-team/PyBaMM/pull/1929)))
+-   Isothermal models now calculate heat source terms (but the temperature remains constant). The models now also account for current collector heating when `dimensionality=0` ([#1929](https://github.com/pybamm-team/PyBaMM/pull/1929))
 -   Added new models for power control and resistance control ([#1917](https://github.com/pybamm-team/PyBaMM/pull/1917))
 -   Initial concentrations can now be provided as a function of `r` as well as `x` ([#1866](https://github.com/pybamm-team/PyBaMM/pull/1866))
 
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index b262d1f4b5..d9be896f26 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -30,6 +30,11 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 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.
+            * "calculate heat source for isothermal models" : str
+                Whether to calculate the heat source terms during isothermal operation.
+                Can be "true" or "false". If "false", the heat source terms are set
+                to zero. Default is "false" since this option may require additional
+                parameters not needed by the electrochemical model.
             * "convection" : str
                 Whether to include the effects of convection in the model. Can be
                 "none" (default), "uniform transverse" or "full transverse".
@@ -169,6 +174,7 @@ def __init__(self, extra_options):
         self.possible_options = {
             "calculate discharge energy": ["false", "true"],
             "cell geometry": ["arbitrary", "pouch"],
+            "calculate heat source for isothermal models": ["false", "true"],
             "convection": ["none", "uniform transverse", "full transverse"],
             "current collector": [
                 "uniform",
diff --git a/pybamm/models/submodels/thermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal.py
index d620c8f38d..2fd8270a71 100644
--- a/pybamm/models/submodels/thermal/isothermal.py
+++ b/pybamm/models/submodels/thermal/isothermal.py
@@ -42,3 +42,39 @@ def get_fundamental_variables(self):
         )
 
         return variables
+
+    def get_coupled_variables(self, variables):
+        if self.options["calculate heat source for isothermal models"] == "true":
+            variables.update(self._get_standard_coupled_variables(variables))
+        else:
+            ieh = "irreversible electrochemical heating"
+            variables.update(
+                {
+                    "Ohmic heating": pybamm.Scalar(0),
+                    "Ohmic heating [W.m-3]": pybamm.Scalar(0),
+                    "X-averaged Ohmic heating": pybamm.Scalar(0),
+                    "X-averaged Ohmic heating [W.m-3]": pybamm.Scalar(0),
+                    "Volume-averaged Ohmic heating": pybamm.Scalar(0),
+                    "Volume-averaged Ohmic heating [W.m-3]": pybamm.Scalar(0),
+                    "Irreversible electrochemical heating": pybamm.Scalar(0),
+                    "Irreversible electrochemical heating [W.m-3]": pybamm.Scalar(0),
+                    "X-averaged " + ieh: pybamm.Scalar(0),
+                    "X-averaged " + ieh + " [W.m-3]": pybamm.Scalar(0),
+                    "Volume-averaged " + ieh: pybamm.Scalar(0),
+                    "Volume-averaged " + ieh + "[W.m-3]": pybamm.Scalar(0),
+                    "Reversible heating": pybamm.Scalar(0),
+                    "Reversible heating [W.m-3]": pybamm.Scalar(0),
+                    "X-averaged reversible heating": pybamm.Scalar(0),
+                    "X-averaged reversible heating [W.m-3]": pybamm.Scalar(0),
+                    "Volume-averaged reversible heating": pybamm.Scalar(0),
+                    "Volume-averaged reversible heating [W.m-3]": pybamm.Scalar(0),
+                    "Total heating": pybamm.Scalar(0),
+                    "Total heating [W.m-3]": pybamm.Scalar(0),
+                    "X-averaged total heating": pybamm.Scalar(0),
+                    "X-averaged total heating [W.m-3]": pybamm.Scalar(0),
+                    "Volume-averaged total heating": pybamm.Scalar(0),
+                    "Volume-averaged total heating [W.m-3]": pybamm.Scalar(0),
+                }
+            )
+
+        return variables
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 542a838498..2a16da3981 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
@@ -16,6 +16,7 @@
 PRINT_OPTIONS_OUTPUT = """\
 'calculate discharge energy': 'false' (possible: ['false', 'true'])
 'cell geometry': 'pouch' (possible: ['arbitrary', 'pouch'])
+'calculate heat source for isothermal models': 'false' (possible: ['false', 'true'])
 'convection': 'none' (possible: ['none', 'uniform transverse', 'full transverse'])
 'current collector': 'uniform' (possible: ['uniform', 'potential pair', 'potential pair quite conductive'])
 'dimensionality': 0 (possible: [0, 1, 2])
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 0fe33e09b7..b19dcbf806 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
@@ -13,6 +13,13 @@ def test_well_posed(self):
         options = {"thermal": "isothermal"}
         self.check_well_posedness(options)
 
+    def test_well_posed_isothermal_heat_source(self):
+        options = {
+            "calculate heat source for isothermal models": "true",
+            "thermal": "isothermal",
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_2plus1D(self):
         options = {"current collector": "potential pair", "dimensionality": 1}
         self.check_well_posedness(options)