diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 18e9b3d1d..b45de016b 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1239,6 +1239,12 @@ def initialise(self):
     def define_temperature(self):
         super().define_temperature()
 
+        # NOTE this won't be needed anymore when https://github.com/FEniCS/dolfinx/pull/4140
+        # is released
+
+        # because dolfinx.fem.Expressions cannot work with submeshes
+        # (ie. mixing parent and submesh),
+        # we need to create "sub" temperature functions for each subdomain
         # pass temperature function to each subdomain
         if isinstance(self.temperature_fenics, fem.Function):
             for subdomain in self.volume_subdomains:
@@ -1866,6 +1872,19 @@ def post_processing(self):
                 export.data.append(c)
                 export.t.append(float(self.t))
 
+    def update_time_dependent_values(self):
+        super().update_time_dependent_values()
+
+        # update sub_T if temperature is given as a function
+        if self.temperature_time_dependent:
+            if isinstance(self.temperature_fenics, fem.Function):
+                for subdomain in self.volume_subdomains:
+                    temp = self.temperature_fenics
+                    sub_T = subdomain.sub_T
+                    from festim.helpers import nmm_interpolate
+
+                    nmm_interpolate(f_out=sub_T, f_in=temp)
+
     def iterate(self):
         """Iterates the model for a given time step."""
         if self.show_progress_bar:
diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py
index bf2326a61..8d678f43f 100644
--- a/test/system_tests/test_misc.py
+++ b/test/system_tests/test_misc.py
@@ -1,3 +1,5 @@
+from mpi4py import MPI
+
 import dolfinx
 import numpy as np
 import pytest
@@ -352,3 +354,38 @@ def test_timesteps():
 
     expected_timesteps = np.linspace(0, 10, num=10, endpoint=False)
     assert np.allclose(my_model.timesteps, expected_timesteps)
+
+
+def test_sub_temperature_as_function_mixed_domain_not_updated():
+
+    my_model = F.HydrogenTransportProblemDiscontinuous()
+
+    n = 8
+    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n)
+    my_model.mesh = F.Mesh(mesh)
+
+    surface_1 = F.SurfaceSubdomain(id=1)
+
+    material_1 = F.Material(D_0=1, E_D=0, name="material_1")
+    volume_1 = F.VolumeSubdomain(id=2, material=material_1)
+
+    my_model.subdomains = [surface_1, volume_1]
+
+    H = F.Species(name="H", subdomains=my_model.volume_subdomains)
+    my_model.species = [H]
+
+    my_model.boundary_conditions = [
+        F.FixedConcentrationBC(species=H, subdomain=surface_1, value=lambda T: 2 * T)
+    ]
+
+    my_model.temperature = lambda t, x: 1 + x[0] + x[1] + t
+
+    my_model.settings = F.Settings(final_time=1, atol=1e-9, rtol=1e-9, stepsize=0.1)
+
+    avg_surf = F.AverageSurface(field=H, surface=surface_1)
+    my_model.exports = [avg_surf]
+
+    my_model.initialise()
+    my_model.run()
+
+    assert not np.allclose(avg_surf.data, 4.0)