Repository: pybamm-team/PyBaMM

Issue:
Add "Surface temperature" submodel
This model accounts for the fact that the surface temperature of a cell in a thermal chamber is typically different from the "set"/"ambient" temperature of the chamber.

Suggested steps to do this below, but open to discussion:
- Add a submodel `SurfaceAmbientTemperature` which sets surface temperature equal to ambient temperature
- Add another submodel `SurfaceTemperatureODE` which implements the equation from [Mohan et al](https://d1wqtxts1xzle7.cloudfront.net/69868813/j.jpowsour.2014.01.09720210918-12240-190ns8h-libre.pdf?1632013230=&response-content-disposition=inline%3B+filename%3DA_lumped_parameter_electro_thermal_model.pdf&Expires=1713470219&Signature=VVdyggn0EMZaZUDAhblutOENZ6faE0Wua0uoq71ebx912paMfn2O2-lJq8SUK1nwkZf7v5jF65dIUF03R~IKZo5LXZdYEOb48oPshUycdG0h3K4EHPaGW07CNjd6e~UFP2JyNjUDMAnfP5oZHD~8cr5MRoTbNi8WTfcJ8SSdc1nyOr4uKy5xevjQKW3D0cvAbu9gSnbVZC52qdwM6jeobcJLAx7AhMetR3x1r2lHd6Jrrz-vznm1JeczXl6EzEzMgqr2tLOcKO7GLhlUVYIDhdGawR76FJHLBNyhAEPXtEMq7IUibswa~~2fdgZNTemj9lCjsExrmrri-xZhHFL5vQ__&Key-Pair-Id=APKAJLOHF5GGSLRBV4ZA):
<img width="347" alt="Screenshot 2024-04-18 at 2 59 00 PM" src="https://github.com/pybamm-team/PyBaMM/assets/20817509/db0728a2-c986-46ff-800e-c5a38b294648">

- In that submodel, register a citation for that paper
- In all other thermal submodels, replace "ambient temperature" with "surface temperature"
- In `BaseBatteryModel`, add an option for "surface temperature" which can be "ambient" (default) or "ODE", then add an if statement for selecting the appropriate model
- Add a test for selecting the "ODE" option
- Add docs

Discussion:
valentinsulzer: Still need to add tests and update examples 
TomTranter: I am happy to review this when ready

valentinsulzer: @TomTranter ready for review
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4203?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
Attention: Patch coverage is `98.50746%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 99.57%. Comparing base [(`088ca38`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/088ca3827c831e1b93c4667b6f4353f64fb32cbf?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`b8a62b2`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/b8a62b2e8a6179fb3764ae138b442af5e23df8ed?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

| [Files](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4203?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) | Patch % | Lines |
|---|---|---|
| [...m/models/full\_battery\_models/base\_battery\_model.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4203?src=pr&el=tree&filepath=pybamm%2Fmodels%2Ffull_battery_models%2Fbase_battery_model.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-cHliYW1tL21vZGVscy9mdWxsX2JhdHRlcnlfbW9kZWxzL2Jhc2VfYmF0dGVyeV9tb2RlbC5weQ==) | 88.88% | [1 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4203?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #4203      +/-   ##
===========================================
- Coverage    99.58%   99.57%   -0.01%     
===========================================
  Files          288      291       +3     
  Lines        21799    21843      +44     
===========================================
+ Hits         21708    21751      +43     
- Misses          91       92       +1     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4203?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

kratman: Looks like the examples are failing here:
```
            for phase in self.options.phases[domain]:
                if particle == "Fickian diffusion":
                    submod = pybamm.particle.FickianDiffusion(
                        self.param, domain, self.options, phase=phase, x_average=True
                    )
                elif particle in [
                    "uniform profile",
                    "quadratic profile",
                    "quartic profile",
                ]:
                    submod = pybamm.particle.XAveragedPolynomialProfile(
                        self.param, domain, self.options, phase=phase
                    )
                elif particle == "MSMR":
                    submod = pybamm.particle.MSMRDiffusion(
                        self.param, domain, self.options, phase=phase, x_average=True
                    )
                    # also set the submodel for calculating stoichiometry from
                    # potential
                    self.submodels[f"{domain} {phase} stoichiometry"] = (
                        pybamm.particle.MSMRStoichiometryVariables(
                            self.param,
                            domain,
                            self.options,
                            phase=phase,
                            x_average=True,
                        )
                    )
                self.submodels[f"{domain} {phase} particle"] = submod
```

The variable submod does not get defined before assignment. PyCharm identifies that the variable may be undefined. We should probably make sure our inspection tools flag this
TomTranter: The model seems to be working but I am not sure about some of the numbers or what the constants in the equation represent. It would make more sense to me to have 3 nodes representing the core, casing and environment. I also would prefer the equations to be consistent with the other heat capacities and specify the thermal conductivity and specific heat capacity. The paper quotes these values for the parameters in the model 
![image](https://github.com/pybamm-team/PyBaMM/assets/7068741/f632f1ba-2731-4b17-a335-ae44bae19988)
The thermal resistance seems ok but capacity seems a bit high unless my assumed casing thickness is a bit low.
```
# Assuming a 21700 cylindrical cell with Al casing

rho_Al = 2710.00  # [kg/m^3]
h_casing = 70.0e-3  # [m] height of casing
t_casing = 1.0e-3  # [m] casing thickness
r_casing = 10.50e-3  # [m] casing radius

# Volume of the cylindrical side
vol_casing_side = np.pi * ((r_casing + t_casing)**2 - r_casing**2) * h_casing  # [m^3]

# Volume of the end caps
vol_casing_caps = 2 * np.pi * r_casing**2 * t_casing  # [m^3]

# Total volume of the casing
vol_casing = vol_casing_side + vol_casing_caps

Cp_Al = 890.0  # [J/kg.K]
C_Al = rho_Al * vol_casing * Cp_Al

print("Casing heat capacity [J/K]:", C_Al)


k_Al = 237 # W.m-1K-1
R_casing = 1 / (k_Al * t_casing)

print("Environment thermal resistance [K/W]:", R_casing)

parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.update(
    {
        "Casing heat capacity [J.K-1]": C_Al,
        "Environment thermal resistance [K.W-1]": R_casing,
    },
    check_already_exists=False,
)
```
Gives
Casing heat capacity [J/K]: 13.33
Environment thermal resistance [K/W]: 4.21

Could we see what adding another equation does to the model to go from casing to cooling channel centre through air. It seems like the fitted parameters for this model might be doing some averaging of properties which is not that obvious to the user
valentinsulzer: I'd like to stick as closely as possible to what is implemented in the referred paper, so that we can say "this is the model from this paper" rather than having to justify choices of equations and parameters. Especially as this model already gives good fit to some data I've been looking at 

Patch:
diff --git a/examples/scripts/compare_surface_temperature.py b/examples/scripts/compare_surface_temperature.py
new file mode 100644
index 000000000..24ccfc5df
--- /dev/null
+++ b/examples/scripts/compare_surface_temperature.py
@@ -0,0 +1,61 @@
+#
+# Compare lithium-ion battery models
+#
+import pybamm
+
+pybamm.set_logging_level("INFO")
+
+# load models
+models = [
+    pybamm.lithium_ion.SPMe(
+        {"thermal": "lumped", "surface temperature": "ambient"},
+        name="ambient surface temperature",
+    ),
+    pybamm.lithium_ion.SPMe(
+        {"thermal": "lumped", "surface temperature": "lumped"},
+        name="lumped surface temperature",
+    ),
+]
+
+experiment = pybamm.Experiment(
+    [
+        "Discharge at 1C until 2.5V",
+        "Rest for 1 hour",
+    ]
+)
+
+parameter_values = pybamm.ParameterValues("Chen2020")
+parameter_values.update(
+    {
+        "Casing heat capacity [J.K-1]": 30,
+        "Environment thermal resistance [K.W-1]": 10,
+    },
+    check_already_exists=False,
+)
+
+# create and run simulations
+sols = []
+for model in models:
+    model.variables["Bulk temperature [°C]"] = (
+        model.variables["Volume-averaged cell temperature [K]"] - 273.15
+    )
+    model.variables["Surface temperature [°C]"] = (
+        model.variables["Surface temperature [K]"] - 273.15
+    )
+    sim = pybamm.Simulation(
+        model, parameter_values=parameter_values, experiment=experiment
+    )
+    sol = sim.solve([0, 3600])
+    sols.append(sol)
+
+# plot
+pybamm.dynamic_plot(
+    sols,
+    [
+        "Voltage [V]",
+        "Bulk temperature [°C]",
+        "Surface temperature [°C]",
+        "Surface total cooling [W]",
+        "Environment total cooling [W]",
+    ],
+)
diff --git a/src/pybamm/CITATIONS.bib b/src/pybamm/CITATIONS.bib
index 8fb9c6dc9..3d853738b 100644
--- a/src/pybamm/CITATIONS.bib
+++ b/src/pybamm/CITATIONS.bib
@@ -252,6 +252,16 @@
   doi = {10.3390/batteries5040064},
 }
 
+@article{lin2014lumped,
+  title={A lumped-parameter electro-thermal model for cylindrical batteries},
+  author={Lin, Xinfan and Perez, Hector E and Mohan, Shankar and Siegel, Jason B and Stefanopoulou, Anna G and Ding, Yi and Castanier, Matthew P},
+  journal={Journal of Power Sources},
+  volume={257},
+  pages={1--11},
+  year={2014},
+  publisher={Elsevier}
+}
+
 @article{Marquis2019,
   title = {{An asymptotic derivation of a single particle model with electrolyte}},
   author = {Marquis, Scott G. and Sulzer, Valentin and Timms, Robert and Please, Colin P. and Chapman, S. Jon},
diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py
index 8a2e44333..ccda594b1 100644
--- a/src/pybamm/models/full_battery_models/base_battery_model.py
+++ b/src/pybamm/models/full_battery_models/base_battery_model.py
@@ -194,6 +194,11 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "surface form" : str
                 Whether to use the surface formulation of the problem. Can be "false"
                 (default), "differential" or "algebraic".
+            * "surface temperature" : str
+                Sets the surface temperature model to use. Can be "ambient" (default),
+                which sets the surface temperature equal to the ambient temperature, or
+                "lumped", which adds an ODE for the surface temperature (e.g. to model
+                internal heating of a thermal chamber).
             * "thermal" : str
                 Sets the thermal model to use. Can be "isothermal" (default), "lumped",
                 "x-lumped", or "x-full". The 'cell geometry' option must be set to
@@ -278,7 +283,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             ],
             "particle": [
                 "Fickian diffusion",
-                "fast diffusion",
                 "uniform profile",
                 "quadratic profile",
                 "quartic profile",
@@ -304,6 +308,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             "SEI porosity change": ["false", "true"],
             "stress-induced diffusion": ["false", "true"],
             "surface form": ["false", "differential", "algebraic"],
+            "surface temperature": ["ambient", "lumped"],
             "thermal": ["isothermal", "lumped", "x-lumped", "x-full"],
             "total interfacial current density as a state": ["false", "true"],
             "transport efficiency": [
@@ -554,16 +559,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 )
 
         # Renamed options
-        if options["particle"] == "fast diffusion":
-            raise pybamm.OptionError(
-                "The 'fast diffusion' option has been renamed. "
-                "Use 'uniform profile' instead."
-            )
-        if options["SEI porosity change"] in [True, False]:
-            raise pybamm.OptionError(
-                "SEI porosity change must now be given in string format "
-                "('true' or 'false')"
-            )
         if options["working electrode"] == "negative":
             raise pybamm.OptionError(
                 "The 'negative' working electrode option has been removed because "
@@ -633,6 +628,12 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                     "be 'none': 'particle mechanics', 'loss of active material'"
                 )
 
+        if options["surface temperature"] == "lumped":
+            if options["thermal"] not in ["isothermal", "lumped"]:
+                raise pybamm.OptionError(
+                    "lumped surface temperature model only compatible with isothermal "
+                    "or lumped thermal model"
+                )
         if "true" in options["SEI on cracks"]:
             sei_on_cr = options["SEI on cracks"]
             p_mechanics = options["particle mechanics"]
@@ -1265,6 +1266,13 @@ class BaseBatteryModel(pybamm.BaseModel):
             self.param, self.options, x_average
         )
 
+    def set_surface_temperature_submodel(self):
+        if self.options["surface temperature"] == "ambient":
+            submodel = pybamm.thermal.surface.Ambient
+        elif self.options["surface temperature"] == "lumped":
+            submodel = pybamm.thermal.surface.Lumped
+        self.submodels["surface temperature"] = submodel(self.param, self.options)
+
     def set_current_collector_submodel(self):
         if self.options["current collector"] in ["uniform"]:
             submodel = pybamm.current_collector.Uniform(self.param)
diff --git a/src/pybamm/models/full_battery_models/lead_acid/full.py b/src/pybamm/models/full_battery_models/lead_acid/full.py
index b5b561e6d..f873c779b 100644
--- a/src/pybamm/models/full_battery_models/lead_acid/full.py
+++ b/src/pybamm/models/full_battery_models/lead_acid/full.py
@@ -27,6 +27,7 @@ class Full(BaseModel):
         self.set_electrolyte_submodel()
         self.set_solid_submodel()
         self.set_thermal_submodel()
+        self.set_surface_temperature_submodel()
         self.set_side_reaction_submodels()
         self.set_current_collector_submodel()
         self.set_sei_submodel()
diff --git a/src/pybamm/models/full_battery_models/lead_acid/loqs.py b/src/pybamm/models/full_battery_models/lead_acid/loqs.py
index c63c9cd11..76a4142a0 100644
--- a/src/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/src/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -27,6 +27,7 @@ class LOQS(BaseModel):
         self.set_electrolyte_submodel()
         self.set_electrode_submodels()
         self.set_thermal_submodel()
+        self.set_surface_temperature_submodel()
         self.set_side_reaction_submodels()
         self.set_current_collector_submodel()
         self.set_sei_submodel()
diff --git a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index 6db56b74c..dfe2512f6 100644
--- a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -48,6 +48,7 @@ class BaseModel(pybamm.BaseBatteryModel):
         self.set_electrolyte_concentration_submodel()
         self.set_electrolyte_potential_submodel()
         self.set_thermal_submodel()
+        self.set_surface_temperature_submodel()
         self.set_current_collector_submodel()
         self.set_sei_submodel()
         self.set_sei_on_cracks_submodel()
diff --git a/src/pybamm/models/submodels/thermal/__init__.py b/src/pybamm/models/submodels/thermal/__init__.py
index cc8f769f3..90bd14bc3 100644
--- a/src/pybamm/models/submodels/thermal/__init__.py
+++ b/src/pybamm/models/submodels/thermal/__init__.py
@@ -2,5 +2,6 @@ from .base_thermal import BaseThermal
 from .isothermal import Isothermal
 from .lumped import Lumped
 from . import pouch_cell
+from . import surface
 
-__all__ = ['base_thermal', 'isothermal', 'lumped', 'pouch_cell']
+__all__ = ["base_thermal", "isothermal", "lumped", "pouch_cell"]
diff --git a/src/pybamm/models/submodels/thermal/isothermal.py b/src/pybamm/models/submodels/thermal/isothermal.py
index edcd47bbd..4b729f129 100644
--- a/src/pybamm/models/submodels/thermal/isothermal.py
+++ b/src/pybamm/models/submodels/thermal/isothermal.py
@@ -73,8 +73,8 @@ class Isothermal(BaseThermal):
                 "Total heating [W]",
                 "Negative current collector Ohmic heating [W.m-3]",
                 "Positive current collector Ohmic heating [W.m-3]",
-                "Lumped total cooling [W.m-3]",
-                "Lumped total cooling [W]",
+                "Surface total cooling [W.m-3]",
+                "Surface total cooling [W]",
             ]:
                 # All variables are zero
                 variables.update({var: zero})
diff --git a/src/pybamm/models/submodels/thermal/lumped.py b/src/pybamm/models/submodels/thermal/lumped.py
index 76af2904b..6915a3b18 100644
--- a/src/pybamm/models/submodels/thermal/lumped.py
+++ b/src/pybamm/models/submodels/thermal/lumped.py
@@ -49,9 +49,9 @@ class Lumped(BaseThermal):
 
         # Newton cooling, accounting for surface area to volume ratio
         T_vol_av = variables["Volume-averaged cell temperature [K]"]
-        T_amb = variables["Volume-averaged ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         V = variables["Cell thermal volume [m3]"]
-        Q_cool_W = -self.param.h_total * (T_vol_av - T_amb) * self.param.A_cooling
+        Q_cool_W = -self.param.h_total * (T_vol_av - T_surf) * self.param.A_cooling
         Q_cool_vol_av = Q_cool_W / V
 
         # Contact resistance heating Q_cr
@@ -67,8 +67,8 @@ class Lumped(BaseThermal):
         variables.update(
             {
                 # Lumped cooling
-                "Lumped total cooling [W.m-3]": Q_cool_vol_av,
-                "Lumped total cooling [W]": Q_cool_W,
+                "Surface total cooling [W.m-3]": Q_cool_vol_av,
+                "Surface total cooling [W]": Q_cool_W,
                 # Contact resistance
                 "Lumped contact resistance heating [W.m-3]": Q_cr_vol_av,
                 "Lumped contact resistance heating [W]": Q_cr_W,
@@ -79,7 +79,7 @@ class Lumped(BaseThermal):
     def set_rhs(self, variables):
         T_vol_av = variables["Volume-averaged cell temperature [K]"]
         Q_vol_av = variables["Volume-averaged total heating [W.m-3]"]
-        Q_cool_vol_av = variables["Lumped total cooling [W.m-3]"]
+        Q_cool_vol_av = variables["Surface total cooling [W.m-3]"]
         Q_cr_vol_av = variables["Lumped contact resistance heating [W.m-3]"]
         rho_c_p_eff_av = variables[
             "Volume-averaged effective heat capacity [J.K-1.m-3]"
diff --git a/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py b/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
index a4908c6f5..fb026a9a0 100644
--- a/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
+++ b/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
@@ -54,7 +54,7 @@ class CurrentCollector1D(BaseThermal):
     def set_rhs(self, variables):
         T_av = variables["X-averaged cell temperature [K]"]
         Q_av = variables["X-averaged total heating [W.m-3]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
@@ -64,13 +64,13 @@ class CurrentCollector1D(BaseThermal):
         cell_volume = self.param.L * self.param.L_y * self.param.L_z
         Q_yz_surface = (
             -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z))
-            * (T_av - T_amb)
+            * (T_av - T_surf)
             * yz_surface_area
             / cell_volume
         )
         Q_edge = (
             -(self.param.h_edge(0, z) + self.param.h_edge(self.param.L_y, z))
-            * (T_av - T_amb)
+            * (T_av - T_surf)
             * edge_area
             / cell_volume
         )
@@ -87,7 +87,7 @@ class CurrentCollector1D(BaseThermal):
 
     def set_boundary_conditions(self, variables):
         param = self.param
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         T_av = variables["X-averaged cell temperature [K]"]
 
         # Find tab locations (top vs bottom)
@@ -118,10 +118,10 @@ class CurrentCollector1D(BaseThermal):
         # Calculate heat fluxes weighted by area
         # Note: can't do y-average of h_edge here since y isn't meshed. Evaluate at
         # midpoint.
-        q_tab_n = -param.n.h_tab * (T_av - T_amb)
-        q_tab_p = -param.p.h_tab * (T_av - T_amb)
-        q_edge_top = -param.h_edge(L_y / 2, L_z) * (T_av - T_amb)
-        q_edge_bottom = -param.h_edge(L_y / 2, 0) * (T_av - T_amb)
+        q_tab_n = -param.n.h_tab * (T_av - T_surf)
+        q_tab_p = -param.p.h_tab * (T_av - T_surf)
+        q_edge_top = -param.h_edge(L_y / 2, L_z) * (T_av - T_surf)
+        q_edge_bottom = -param.h_edge(L_y / 2, 0) * (T_av - T_surf)
         q_top = (
             q_tab_n * neg_tab_area * neg_tab_top_bool
             + q_tab_p * pos_tab_area * pos_tab_top_bool
diff --git a/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py b/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
index 7955ee4c3..b2d69ff1b 100644
--- a/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
+++ b/src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
@@ -54,15 +54,15 @@ class CurrentCollector2D(BaseThermal):
     def set_rhs(self, variables):
         T_av = variables["X-averaged cell temperature [K]"]
         Q_av = variables["X-averaged total heating [W.m-3]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
         # Calculate cooling
         Q_yz_surface_W_per_m2 = -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z)) * (
-            T_av - T_amb
+            T_av - T_surf
         )
-        Q_edge_W_per_m2 = -self.param.h_edge(y, z) * (T_av - T_amb)
+        Q_edge_W_per_m2 = -self.param.h_edge(y, z) * (T_av - T_surf)
 
         # Account for surface area to volume ratio of pouch cell in surface cooling
         # term
@@ -98,14 +98,14 @@ class CurrentCollector2D(BaseThermal):
 
     def set_boundary_conditions(self, variables):
         T_av = variables["X-averaged cell temperature [K]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
         # Calculate heat fluxes
-        q_tab_n = -self.param.n.h_tab * (T_av - T_amb)
-        q_tab_p = -self.param.p.h_tab * (T_av - T_amb)
-        q_edge = -self.param.h_edge(y, z) * (T_av - T_amb)
+        q_tab_n = -self.param.n.h_tab * (T_av - T_surf)
+        q_tab_p = -self.param.p.h_tab * (T_av - T_surf)
+        q_edge = -self.param.h_edge(y, z) * (T_av - T_surf)
 
         # Subtract the edge cooling from the tab portion so as to not double count
         # Note: tab cooling is also only applied on the current collector hence
diff --git a/src/pybamm/models/submodels/thermal/pouch_cell/x_full.py b/src/pybamm/models/submodels/thermal/pouch_cell/x_full.py
index f4aa07c56..630ce94b0 100644
--- a/src/pybamm/models/submodels/thermal/pouch_cell/x_full.py
+++ b/src/pybamm/models/submodels/thermal/pouch_cell/x_full.py
@@ -79,7 +79,7 @@ class OneDimensionalX(BaseThermal):
         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]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
@@ -140,28 +140,28 @@ class OneDimensionalX(BaseThermal):
                 (
                     pybamm.boundary_value(lambda_n, "left")
                     * pybamm.boundary_gradient(T_n, "left")
-                    - h_cn * (T_cn - T_amb)
+                    - h_cn * (T_cn - T_surf)
                 )
                 / L_cn
                 + Q_cn
-                + cooling_coefficient_cn * (T_cn - T_amb)
+                + cooling_coefficient_cn * (T_cn - T_surf)
             )
             / self.param.n.rho_c_p_cc(T_cn),
             T: (
                 pybamm.div(lambda_ * pybamm.grad(T))
                 + Q
-                + cooling_coefficient * (T - T_amb)
+                + cooling_coefficient * (T - T_surf)
             )
             / rho_c_p,
             T_cp: (
                 (
                     -pybamm.boundary_value(lambda_p, "right")
                     * pybamm.boundary_gradient(T_p, "right")
-                    - h_cp * (T_cp - T_amb)
+                    - h_cp * (T_cp - T_surf)
                 )
                 / L_cp
                 + Q_cp
-                + cooling_coefficient_cp * (T_cp - T_amb)
+                + cooling_coefficient_cp * (T_cp - T_surf)
             )
             / self.param.p.rho_c_p_cc(T_cp),
         }
diff --git a/src/pybamm/models/submodels/thermal/surface/__init__.py b/src/pybamm/models/submodels/thermal/surface/__init__.py
new file mode 100644
index 000000000..dece4b44a
--- /dev/null
+++ b/src/pybamm/models/submodels/thermal/surface/__init__.py
@@ -0,0 +1,2 @@
+from .ambient import Ambient
+from .lumped import Lumped
diff --git a/src/pybamm/models/submodels/thermal/surface/ambient.py b/src/pybamm/models/submodels/thermal/surface/ambient.py
new file mode 100644
index 000000000..a28591af3
--- /dev/null
+++ b/src/pybamm/models/submodels/thermal/surface/ambient.py
@@ -0,0 +1,33 @@
+#
+# Class for ambient surface temperature submodel
+import pybamm
+
+
+class Ambient(pybamm.BaseSubModel):
+    """
+    Class for setting surface temperature equal to ambient temperature.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param, options=options)
+
+    def get_coupled_variables(self, variables):
+        T_amb = variables["Ambient temperature [K]"]
+        T_amb_av = variables["Volume-averaged ambient temperature [K]"]
+
+        variables.update(
+            {
+                "Surface temperature [K]": T_amb,
+                "Volume-averaged surface temperature [K]": T_amb_av,
+                "Environment total cooling [W]": pybamm.Scalar(0),
+            }
+        )
+        return variables
diff --git a/src/pybamm/models/submodels/thermal/surface/lumped.py b/src/pybamm/models/submodels/thermal/surface/lumped.py
new file mode 100644
index 000000000..dc481947e
--- /dev/null
+++ b/src/pybamm/models/submodels/thermal/surface/lumped.py
@@ -0,0 +1,51 @@
+#
+# Class for ambient surface temperature submodel
+import pybamm
+
+
+class Lumped(pybamm.BaseSubModel):
+    """
+    Class for the lumped surface temperature submodel, which adds an ODE for the
+    surface temperature.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param, options=options)
+        pybamm.citations.register("lin2014lumped")
+
+    def get_fundamental_variables(self):
+        T_surf = pybamm.Variable("Surface temperature [K]")
+        variables = {"Surface temperature [K]": T_surf}
+
+        return variables
+
+    def get_coupled_variables(self, variables):
+        T_surf = variables["Surface temperature [K]"]
+        T_amb = variables["Ambient temperature [K]"]
+        R_env = pybamm.Parameter("Environment thermal resistance [K.W-1]")
+        Q_cool_env = -(T_surf - T_amb) / R_env
+        variables["Environment total cooling [W]"] = Q_cool_env
+        return variables
+
+    def set_rhs(self, variables):
+        T_surf = variables["Surface temperature [K]"]
+
+        Q_cool_bulk = variables["Surface total cooling [W]"]
+        Q_heat_bulk = -Q_cool_bulk
+
+        Q_cool_env = variables["Environment total cooling [W]"]
+        rho_c_p_case = pybamm.Parameter("Casing heat capacity [J.K-1]")
+
+        self.rhs[T_surf] = (Q_heat_bulk + Q_cool_env) / rho_c_p_case
+
+    def set_initial_conditions(self, variables):
+        T_surf = variables["Surface temperature [K]"]
+        self.initial_conditions = {T_surf: self.param.T_init}
diff --git a/src/pybamm/plotting/plot_thermal_components.py b/src/pybamm/plotting/plot_thermal_components.py
index 8cb99a454..e45a08112 100644
--- a/src/pybamm/plotting/plot_thermal_components.py
+++ b/src/pybamm/plotting/plot_thermal_components.py
@@ -55,7 +55,7 @@ def plot_thermal_components(
     volume = solution["Cell thermal volume [m3]"].entries
 
     heating_sources = [
-        "Lumped total cooling",
+        "Surface total cooling",
         "Ohmic heating",
         "Irreversible electrochemical heating",
         "Reversible heating",
@@ -77,9 +77,9 @@ def plot_thermal_components(
     # Plot
     # Initialise
     total_heat = 0
-    bottom_heat = heats["Lumped total cooling"]
+    bottom_heat = heats["Surface total cooling"]
     total_cumul_heat = 0
-    bottom_cumul_heat = cumul_heats["Lumped total cooling"]
+    bottom_cumul_heat = cumul_heats["Surface total cooling"]
     # Plot components
     for name in heating_sources:
         top_heat = bottom_heat + abs(heats[name])


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
index 603f64d71..ac58f40bb 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
@@ -101,6 +101,35 @@ class TestThermal:
 
         assert 1e-5 > err(solutions["SPMe 1+1D"], solutions["SPMe 2+1D"])
 
+    def test_surface_temperature_models(self):
+        models = {
+            option: pybamm.lithium_ion.SPM(
+                {"thermal": "lumped", "surface temperature": option}
+            )
+            for option in ["lumped", "ambient"]
+        }
+
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values.update(
+            {
+                "Casing heat capacity [J.K-1]": 30,
+                "Environment thermal resistance [K.W-1]": 10,
+            },
+            check_already_exists=False,
+        )
+
+        sols = {}
+        for name, model in models.items():
+            sim = pybamm.Simulation(model, parameter_values=parameter_values)
+            sol = sim.solve([0, 3600])
+            sols[name] = sol
+
+        for var in ["Volume-averaged cell temperature [K]", "Surface temperature [K]"]:
+            # ignore first entry as it is the initial condition
+            T_ambient_model = sols["ambient"][var].entries[1:]
+            T_lumped_model = sols["lumped"][var].entries[1:]
+            np.testing.assert_array_less(T_ambient_model, T_lumped_model)
+
     def test_lumped_contact_resistance(self):
         # Test that the heating with contact resistance is greater than without
 
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 6ee38faf9..033dcf534 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
@@ -36,7 +36,7 @@ PRINT_OPTIONS_OUTPUT = """\
 'number of MSMR reactions': 'none' (possible: ['none'])
 'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk'])
 'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'differential power', 'explicit power', 'resistance', 'differential resistance', 'explicit resistance', 'CCCV'])
-'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'fast diffusion', 'uniform profile', 'quadratic profile', 'quartic profile', 'MSMR'])
+'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'uniform profile', 'quadratic profile', 'quartic profile', 'MSMR'])
 'particle mechanics': 'swelling only' (possible: ['none', 'swelling only', 'swelling and cracking'])
 'particle phases': '1' (possible: ['1', '2'])
 'particle shape': 'spherical' (possible: ['spherical', 'no particles'])
@@ -47,6 +47,7 @@ PRINT_OPTIONS_OUTPUT = """\
 'SEI porosity change': 'false' (possible: ['false', 'true'])
 'stress-induced diffusion': 'true' (possible: ['false', 'true'])
 'surface form': 'differential' (possible: ['false', 'differential', 'algebraic'])
+'surface temperature': 'ambient' (possible: ['ambient', 'lumped'])
 'thermal': 'x-full' (possible: ['isothermal', 'lumped', 'x-lumped', 'x-full'])
 'total interfacial current density as a state': 'false' (possible: ['false', 'true'])
 'transport efficiency': 'Bruggeman' (possible: ['Bruggeman', 'ordered packing', 'hyperbola of revolution', 'overlapping spheres', 'tortuosity factor', 'random overlapping cylinders', 'heterogeneous catalyst', 'cation-exchange membrane'])
@@ -211,8 +212,6 @@ class TestBaseBatteryModel(unittest.TestCase):
             pybamm.BaseBatteryModel({"convection": "full transverse"})
         with self.assertRaisesRegex(pybamm.OptionError, "particle"):
             pybamm.BaseBatteryModel({"particle": "bad particle"})
-        with self.assertRaisesRegex(pybamm.OptionError, "The 'fast diffusion'"):
-            pybamm.BaseBatteryModel({"particle": "fast diffusion"})
         with self.assertRaisesRegex(pybamm.OptionError, "working electrode"):
             pybamm.BaseBatteryModel({"working electrode": "bad working electrode"})
         with self.assertRaisesRegex(pybamm.OptionError, "The 'negative' working"):
@@ -233,10 +232,6 @@ class TestBaseBatteryModel(unittest.TestCase):
             pybamm.BaseBatteryModel({"SEI film resistance": "bad SEI film resistance"})
         with self.assertRaisesRegex(pybamm.OptionError, "SEI porosity change"):
             pybamm.BaseBatteryModel({"SEI porosity change": "bad SEI porosity change"})
-        with self.assertRaisesRegex(
-            pybamm.OptionError, "SEI porosity change must now be given in string format"
-        ):
-            pybamm.BaseBatteryModel({"SEI porosity change": True})
         # changing defaults based on other options
         model = pybamm.BaseBatteryModel()
         self.assertEqual(model.options["SEI film resistance"], "none")
@@ -387,6 +382,12 @@ class TestBaseBatteryModel(unittest.TestCase):
                 }
             )
 
+        # surface thermal model
+        with self.assertRaisesRegex(pybamm.OptionError, "surface temperature"):
+            pybamm.BaseBatteryModel(
+                {"surface temperature": "lumped", "thermal": "x-full"}
+            )
+
         # phases
         with self.assertRaisesRegex(pybamm.OptionError, "multiple particle phases"):
             pybamm.BaseBatteryModel({"particle phases": "2", "surface form": "false"})
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 c8a3f6b50..9c093c0c6 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
@@ -31,6 +31,10 @@ class BaseUnitTestLithiumIon:
         options = {"thermal": "lumped"}
         self.check_well_posedness(options)
 
+    def test_well_posed_lumped_thermal_model_surface_temperature(self):
+        options = {"thermal": "lumped", "surface temperature": "lumped"}
+        self.check_well_posedness(options)
+
     def test_well_posed_x_full_thermal_model(self):
         options = {"thermal": "x-full"}
         self.check_well_posedness(options)