diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index d3f92714d..0c052a22d 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -1173,10 +1173,6 @@ class SNES_Scalar(SolverBaseClass):
 
 
         self.essential_bcs = []
-        # TODO(BUG): add_natural_bc() causes PETSc error 73 ("Object in wrong state")
-        # when used with this solver. The Stokes solver's natural BCs work correctly,
-        # suggesting a setup/ordering issue specific to scalar Poisson.
-        # See planning file: underworld.md (Bugs section, 2026-01-19)
         self.natural_bcs = []
         self.bcs = self.essential_bcs
         self.boundary_conditions = False
@@ -1321,7 +1317,7 @@ class SNES_Scalar(SolverBaseClass):
             bc = PetscDSAddBoundary_UW(cdm.dm,
                                 bc_type,
                                 str(boundary+f"{bc.components}").encode('utf8'),
-                                str(boundary).encode('utf8'),
+                                "UW_Boundaries".encode('utf8'),   # consolidated boundary label
                                 bc.f_id,  # field ID in the DM
                                 num_constrained_components,
                                 <const PetscInt *> &comps_view[0],
@@ -1540,7 +1536,46 @@ class SNES_Scalar(SolverBaseClass):
 
         ## Now add the boundary residual / jacobian terms
 
+        cdef DMLabel c_label
+
+        for bc in self.natural_bcs:
+
+            boundary = bc.boundary
+            boundary_id = bc.PETScID
+
+            value = self.mesh.boundaries[bc.boundary].value
+            bc_label = self.dm.getLabel("UW_Boundaries")
+
+            label_val = value
+
+            i_bd_res = self.ext_dict.bd_res
+            i_bd_jac = self.ext_dict.bd_jac
+
+            c_label = bc_label
 
+            if bc.fn_f is not None:
+
+                UW_PetscDSSetBdResidual(ds.ds, c_label.dmlabel, label_val, boundary_id,
+                                0, 0,
+                                ext.fns_bd_residual[i_bd_res[bc.fns["u_f0"]]],
+                                NULL,
+                                )
+
+                UW_PetscDSSetBdJacobian(ds.ds, c_label.dmlabel, label_val, boundary_id,
+                                0, 0, 0,
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G0"]]],
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G1"]]],
+                                NULL,
+                                NULL,
+                                )
+
+                UW_PetscDSSetBdJacobianPreconditioner(ds.ds, c_label.dmlabel, label_val, boundary_id,
+                                0, 0, 0,
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G0"]]],
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G1"]]],
+                                NULL,
+                                NULL,
+                                )
 
         # Rebuild this lot
 
diff --git a/tests/test_1000_poissonNaturalBC.py b/tests/test_1000_poissonNaturalBC.py
index e14944af9..f3b26cc02 100644
--- a/tests/test_1000_poissonNaturalBC.py
+++ b/tests/test_1000_poissonNaturalBC.py
@@ -53,7 +53,6 @@
 )
 
 
-@pytest.mark.skip(reason="Poisson natural BC setup failing - needs solver investigation (PETSc error 73: Object in wrong state)")
 @pytest.mark.parametrize("mesh", [mesh_simp_reg, mesh_simp_irreg, mesh_quad])
 def test_poisson_natural_bc(mesh):
     """Test Poisson solver with natural (flux) boundary conditions."""