diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index b971d39bb..0cf83440f 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1147,6 +1147,7 @@ def initialise(self):
         self.define_temperature()
         self.convert_source_input_values_to_fenics_objects()
         self.convert_advection_term_to_fenics_objects()
+        self.define_boundary_conditions()
         self.create_flux_values_fenics()
         self.create_initial_conditions()
 
@@ -1154,7 +1155,6 @@ def initialise(self):
             self.create_subdomain_formulation(subdomain)
             subdomain.u.name = f"u_{subdomain.id}"
 
-        self.define_boundary_conditions()
         self.create_formulation()
         self.create_solver()
         self.initialise_exports()
@@ -1408,7 +1408,7 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
                 )
 
         # add fluxes
-        for bc in self.boundary_conditions:
+        for bc in self._unpacked_bcs:
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
                 # check that the bc is applied on a surface
                 # belonging to this subdomain
@@ -1631,7 +1631,7 @@ def create_solver(self):
 
     def create_flux_values_fenics(self):
         """For each particle flux create the ``value_fenics`` attribute"""
-        for bc in self.boundary_conditions:
+        for bc in self._unpacked_bcs:
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
                 volume_subdomain = self.surface_to_volume[bc.subdomain]
                 bc.create_value_fenics(
@@ -1908,7 +1908,7 @@ def create_formulation(self):
             )
 
         # add fluxes
-        for bc in self.boundary_conditions:
+        for bc in self._unpacked_bcs:
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
                 self.formulation -= (
                     bc.value_fenics
diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py
index 393ccb511..95de032d7 100644
--- a/test/test_h_transport_problem.py
+++ b/test/test_h_transport_problem.py
@@ -1373,3 +1373,87 @@ def test_traps_with_CG_elements():
 
     # TEST
     assert my_model.function_space.sub(0).element.basix_element.discontinuous is False
+
+
+def test_surface_reaction_BC_discontinuous():
+    """Test that surface reaction BC conserves flux"""
+
+    # BUILD
+    my_model = F.HydrogenTransportProblemDiscontinuous()
+
+    right_mat = F.Material(
+        D_0=1,
+        E_D=0,
+        K_S_0=1,
+        E_K_S=0,
+    )
+    left_mat = F.Material(D_0=1, E_D=0, K_S_0=2, E_K_S=0)
+
+    left_vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=left_mat)
+    right_vol = F.VolumeSubdomain1D(
+        id=2,
+        borders=[1, 2],
+        material=right_mat,
+    )
+
+    vertices_left = np.linspace(0, 1, num=1000, endpoint=False)
+    vertices_right = np.linspace(1, 2, num=1000)
+    vertices = np.concatenate([vertices_left, vertices_right])
+    my_model.mesh = F.Mesh1D(vertices)
+
+    inlet = F.SurfaceSubdomain1D(id=3, x=0)
+    outlet = F.SurfaceSubdomain1D(id=4, x=2)
+
+    my_model.subdomains = [
+        left_vol,
+        right_vol,
+        inlet,
+        outlet,
+    ]
+
+    my_model.surface_to_volume = {
+        inlet: left_vol,
+        outlet: right_vol,
+    }
+
+    my_model.method_interface = F.InterfaceMethod.nitsche
+    my_model.interfaces = [
+        F.Interface(id=5, subdomains=[left_vol, right_vol], penalty_term=10),
+    ]
+
+    H = F.Species("H", subdomains=my_model.volume_subdomains)
+    my_model.species = [H]
+
+    my_model.temperature = 600
+
+    my_model.boundary_conditions = [
+        F.FixedConcentrationBC(
+            subdomain=inlet,
+            value=1,
+            species=H,
+        ),
+        F.SurfaceReactionBC(
+            reactant=[H, H],
+            gas_pressure=0,
+            k_r0=1,
+            E_kr=0,
+            k_d0=0,
+            E_kd=0,
+            subdomain=outlet,
+        ),
+    ]
+    permeation_flux = F.SurfaceFlux(field=H, surface=outlet)
+    inlet_flux = F.SurfaceFlux(field=H, surface=inlet)
+    my_model.exports = [permeation_flux, inlet_flux]
+
+    my_model.settings = F.Settings(
+        atol=1e-10,
+        rtol=1e-09,
+        transient=False,
+    )
+
+    my_model.initialise()
+    my_model.run()
+
+    # TEST
+    assert np.isclose(permeation_flux.data[-1], -inlet_flux.data[-1], rtol=1e-5)