diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index 79ece869..13675c49 100644
--- a/docs/source/expression_tree/unary_operator.rst
+++ b/docs/source/expression_tree/unary_operator.rst
@@ -25,6 +25,11 @@ Unary Operators
 .. autoclass:: pybamm.NumpyBroadcast
   :members:
 
+.. autoclass:: pybamm.SurfaceValue
+  :members:
+
 .. autofunction:: pybamm.grad
 
 .. autofunction:: pybamm.div
+
+.. autofunction:: pybamm.surf
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index ac4b6da7..98c5b202 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -78,8 +78,10 @@ from .expression_tree.unary_operators import (
     Divergence,
     Broadcast,
     NumpyBroadcast,
+    SurfaceValue,
     grad,
     div,
+    surf,
 )
 from .expression_tree.function_parameter import FunctionParameter
 from .expression_tree.scalar import Scalar
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index a572b4fb..18cff23e 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -212,6 +212,13 @@ class Discretisation(object):
                 )
             return symbol
 
+        elif isinstance(symbol, pybamm.SurfaceValue):
+            child = symbol.children[0]
+            discretised_child = self.process_symbol(child)
+            return self._spatial_methods[symbol.domain[0]].surface_value(
+                discretised_child
+            )
+
         elif isinstance(symbol, pybamm.BinaryOperator):
             return self.process_binary_operators(symbol)
 
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 17ce41fd..abe2babb 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -210,6 +210,16 @@ class NumpyBroadcast(Broadcast):
             return child_eval * self.broadcasting_vector
 
 
+class SurfaceValue(SpatialOperator):
+    """A node in the expression tree which gets the surface value of a variable.
+
+    **Extends:** :class:`SpatialOperator`
+    """
+
+    def __init__(self, child):
+        super().__init__("surf", child)
+
+
 #
 # Methods to call Gradient and Divergence
 #
@@ -251,3 +261,27 @@ def div(expression):
     """
 
     return Divergence(expression)
+
+
+#
+# Method to call SurfaceValue
+#
+
+
+def surf(variable):
+    """convenience function for creating a :class:`SurfaceValue`
+
+    Parameters
+    ----------
+
+    variable : :class:`Symbol`
+        the surface value of this variable will be returned
+
+    Returns
+    -------
+
+    :class:`GetSurfaceValue`
+        the surface value of ``variable``
+    """
+
+    return SurfaceValue(variable)
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index b172ad3d..86cc3f68 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -237,6 +237,33 @@ class FiniteVolume(pybamm.SpatialMethod):
             left_ghost_cell, discretised_symbol, right_ghost_cell
         )
 
+    def surface_value(self, discretised_symbol):
+        """
+        Uses linear extrapolation to get the surface value of a variable in the
+        Finite Volume Method.
+
+        Parameters
+        -----------
+        discretised_symbol : :class:`pybamm.StateVector`
+            The discretised variable (a state vector) from which to calculate
+            the surface value.
+
+        Returns
+        -------
+        :class:`pybamm.Variable`
+            The variable representing the surface value.
+        """
+        # Better to make class similar NodeToEdge and pass function?
+        # def surface_value(array):
+        #     "Linear extrapolation for surface value"
+        #     array[-1] + (array[-1] - array[-2]) / 2
+        # ... or make StateVector and add?
+        y_slice_stop = discretised_symbol.y_slice.stop
+        last_node = pybamm.StateVector(slice(y_slice_stop - 1, y_slice_stop))
+        penultimate_node = pybamm.StateVector(slice(y_slice_stop - 2, y_slice_stop - 1))
+        surface_value = (last_node + (last_node - penultimate_node) / 2)
+        return surface_value
+
     #######################################################
     # Can probably be moved outside of the spatial method
     ######################################################
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 1006ed32..b92104b7 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -102,6 +102,24 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
+    def surface_value(self, discretised_symbol):
+        """
+        Returns the surface value using the approriate expression for the
+        spatial method.
+
+        Parameters
+        -----------
+        discretised_symbol : :class:`pybamm.StateVector`
+            The discretised variable (a state vector) from which to calculate
+            the surface value.
+
+        Returns
+        -------
+        :class:`pybamm.Variable`
+            The variable representing the surface value.
+        """
+        raise NotImplementedError
+
     # We could possibly move the following outside of SpatialMethod
     # depending on the requirements of the FiniteVolume
 

Get surface values
**Summary**
Implement a class which takes a variable and returns the value at the surface. 

**Motivation**
This is required for implementing the SPM #72 as the concentrations on the surface of the particles are needed in the Butler-Volmer expression and the OCPs.
diff --git a/tests/test_spatial_methods/test_finite_volume.py b/tests/test_spatial_methods/test_finite_volume.py
index 7c7a2e71..ce84a225 100644
--- a/tests/test_spatial_methods/test_finite_volume.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -32,6 +32,31 @@ class TestFiniteVolume(unittest.TestCase):
         avd = pybamm.NodeToEdge(d, arithmetic_mean)
         np.testing.assert_array_equal(avd.evaluate(None, y_test), np.ones(9))
 
+    def test_surface_value(self):
+        # create discretisation
+        defaults = shared.TestDefaults1DParticle(10)
+        spatial_methods = {"negative particle": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+        mesh = disc.mesh
+
+        combined_submesh = mesh.combine_submeshes("negative particle")
+
+        # create variable
+        var = pybamm.Variable("var", domain="negative particle")
+        surf_eqn = pybamm.surf(var)
+        disc._variables = [var]
+        disc.set_variable_slices()
+        surf_eqn_disc = disc.process_symbol(surf_eqn)
+
+        # check constant extrapolates to constant
+        constant_y = np.ones_like(combined_submesh.nodes)
+        self.assertEqual(surf_eqn_disc.evaluate(None, constant_y), 1)
+
+        # check linear variable extrapolates correctly
+        linear_y = combined_submesh.nodes
+        y_surf = combined_submesh.nodes[-1] + combined_submesh.d_nodes[-1] / 2
+        self.assertEqual(surf_eqn_disc.evaluate(None, linear_y), y_surf)
+
     def test_discretise_diffusivity_times_spatial_operator(self):
         # Set up
         whole_cell = ["negative electrode", "separator", "positive electrode"]

pip install -e .[dev,docs]
parse_log_pytest
pytest
pytest
apt-get update apt-get install -y gcc gfortran libopenblas-dev liblapack-dev
3.6
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs @ file:///opt/conda/conda-bld/attrs_1642510447205/work
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
flake8==5.0.4
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
more-itertools @ file:///tmp/build/80754af9/more-itertools_1637733554872/work
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging @ file:///tmp/build/80754af9/packaging_1637314298585/work
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy @ file:///tmp/build/80754af9/pluggy_1615976315926/work
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py @ file:///opt/conda/conda-bld/py_1644396412707/work
-e git+https://github.com/tinosulzer/PyBaMM.git@b6a3e58f579b1e2b87bb97854e3c902c189e0a7f#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing @ file:///tmp/build/80754af9/pyparsing_1635766073266/work
pyrsistent==0.18.0
pytest==6.2.4
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
toml @ file:///tmp/build/80754af9/toml_1616166611790/work
tornado==6.1
traitlets==4.3.3
typing_extensions @ file:///opt/conda/conda-bld/typing_extensions_1647553014482/work
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp @ file:///tmp/build/80754af9/zipp_1633618647012/work

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - attrs=21.4.0=pyhd3eb1b0_0
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - importlib_metadata=4.8.1=hd3eb1b0_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - more-itertools=8.12.0=pyhd3eb1b0_0
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - packaging=21.3=pyhd3eb1b0_0
  - pip=21.2.2=py36h06a4308_0
  - pluggy=0.13.1=py36h06a4308_0
  - py=1.11.0=pyhd3eb1b0_0
  - pyparsing=3.0.4=pyhd3eb1b0_0
  - pytest=6.2.4=py36h06a4308_2
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - toml=0.10.2=pyhd3eb1b0_0
  - typing_extensions=4.1.1=pyh06a4308_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zipp=3.6.0=pyhd3eb1b0_0
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - flake8==5.0.4
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyrsistent==0.18.0
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tornado==6.1
      - traitlets==4.3.3
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
prefix: /opt/conda/envs/PyBaMM

tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_surface_value
tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_add_ghost_nodes tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_discretise_diffusivity_times_spatial_operator tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_div_convergence tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_div_convergence_internal tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_convergence_with_bcs tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_convergence_without_bcs tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_bcs tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Neumann_bcs tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_mixed_domain tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_node_to_edge tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Dirichlet_bcs tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Neumann_bcs tests/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_operators
pybamm/__init__.py pybamm/expression_tree/unary_operators.py docs/source/expression_tree/unary_operator.rst pybamm/spatial_methods/finite_volume.py pybamm/discretisations/discretisation.py pybamm/spatial_methods/spatial_method.py
pybamm/__init__.py pybamm/expression_tree/unary_operators.py docs/source/expression_tree/unary_operator.rst pybamm/spatial_methods/finite_volume.py pybamm/discretisations/discretisation.py pybamm/spatial_methods/spatial_method.py