diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index 5aaf7ea123..c1a4e5bbfe 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/src/pybamm/spatial_methods/finite_volume.py
@@ -128,6 +128,54 @@ def gradient_matrix(self, domain, domains):
 
         return pybamm.Matrix(matrix)
 
+    def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
+        """
+        Computes the square of the gradient of a symbol.
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.Symbol`
+            The symbol to compute the square of the gradient for.
+        discretised_symbol : :class:`pybamm.Vector`
+            The discretised variable to compute the square of the gradient for.
+        boundary_conditions : dict
+            Boundary conditions for the symbol.
+
+        Returns
+        -------
+        float
+            The gradient squared of the symbol.
+        """
+        domain = symbol.domain
+
+        if symbol in boundary_conditions:
+            bcs = boundary_conditions[symbol]
+            if any(bc[1] == "Dirichlet" for bc in bcs.values()):
+                discretised_symbol, domain = self.add_ghost_nodes(
+                    symbol, discretised_symbol, bcs
+                )
+            elif any(bc[1] == "Neumann" for bc in bcs.values()):
+                # Neumann boundary conditions will be applied after computing the gradient squared
+                pass
+
+        gradient_matrix = self.gradient_matrix(domain, symbol.domains)
+
+        # Compute gradient squared matrix: (∇u)^2 = u^T (L^T L) u
+        gradient_squared_matrix = gradient_matrix.T @ gradient_matrix
+        gradient_squared_result = (
+            discretised_symbol.T @ gradient_squared_matrix @ discretised_symbol
+        )
+
+        # Add Neumann boundary conditions if defined
+        if symbol in boundary_conditions and any(
+            bc[1] == "Neumann" for bc in bcs.values()
+        ):
+            gradient_squared_result = self.add_neumann_values(
+                symbol, gradient_squared_result, bcs, domain
+            )
+
+        return gradient_squared_result.item()
+
     def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the divergence operator.
         See :meth:`pybamm.SpatialMethod.divergence`
diff --git a/tests/unit/test_spatial_methods/test_scikit_finite_element.py b/tests/unit/test_spatial_methods/test_scikit_finite_element.py
index 42b282e08a..2244be4f57 100644
--- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py
+++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py
@@ -157,6 +157,54 @@ def test_gradient(self):
         ans = eqn_disc.evaluate(None, 3 * y**2)
         np.testing.assert_array_less(0, ans)
 
+    def test_gradient_squared(self):
+        mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32, include_particles=False)
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume(),
+            "current collector": pybamm.ScikitFiniteElement(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+
+        y = mesh["current collector"].coordinates[0, :]
+        z = mesh["current collector"].coordinates[1, :]
+
+        gradient = pybamm.grad(var)
+        grad_disc = disc.process_symbol(gradient)
+        grad_disc_y, grad_disc_z = grad_disc.children
+
+        gradient_squared = pybamm.grad_squared(var)
+        gradient_squared_disc = disc.process_symbol(gradient_squared)
+
+        inner_product_y = grad_disc_y.evaluate(None, 5 * y + 6 * z)
+        inner_product_z = grad_disc_z.evaluate(None, 5 * y + 6 * z)
+        inner_product = inner_product_y**2 + inner_product_z**2
+
+        grad_squared_result = gradient_squared_disc.evaluate(None, 5 * y + 6 * z)
+        np.testing.assert_array_almost_equal(grad_squared_result, inner_product)
+        np.testing.assert_array_less(0, grad_squared_result)
+
+        gradient_squared_disc_dirichlet = disc.process_symbol(gradient_squared)
+        grad_squared_result_dirichlet = gradient_squared_disc_dirichlet.evaluate(
+            None, 5 * y + 6 * z
+        )
+
+        np.testing.assert_array_almost_equal(
+            grad_squared_result_dirichlet, inner_product
+        )
+        np.testing.assert_array_less(0, grad_squared_result_dirichlet)
+
+        gradient_squared_disc_neumann = disc.process_symbol(gradient_squared)
+        grad_squared_result_neumann = gradient_squared_disc_neumann.evaluate(
+            None, 5 * y + 6 * z
+        )
+
+        np.testing.assert_array_almost_equal(grad_squared_result_neumann, inner_product)
+
+        np.testing.assert_array_less(0, grad_squared_result_neumann)
+
     def test_manufactured_solution(self):
         mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32, include_particles=False)
         spatial_methods = {