diff --git a/CHANGELOG.md b/CHANGELOG.md
index ef0e0448c6..9425dc6025 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- Allow the option "surface form" to be "differential" in the `MPM` ([#2533](https://github.com/pybamm-team/PyBaMM/pull/2533))
 - Added variables "Loss of lithium due to loss of active material in negative/positive electrode [mol]". These should be included in the calculation of "total lithium in system" to make sure that lithium is truly conserved. ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
 - `initial_soc` can now be a string "x V", in which case the simulation is initialized to start from that voltage ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 - The `ElectrodeSOH` solver can now calculate electrode balance based on a target "cell capacity" (requires cell capacity "Q" as input), as well as the default "cyclable cell capacity" (requires cyclable lithium capacity "Q_Li" as input). Use the keyword argument `known_value` to control which is used. ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
diff --git a/pybamm/models/full_battery_models/lithium_ion/mpm.py b/pybamm/models/full_battery_models/lithium_ion/mpm.py
index 1320d72932..5a171b3f11 100644
--- a/pybamm/models/full_battery_models/lithium_ion/mpm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/mpm.py
@@ -40,41 +40,28 @@ class MPM(SPM):
     """
 
     def __init__(self, options=None, name="Many-Particle Model", build=True):
-        # Necessary options
-        if options is None:
-            options = {"particle size": "distribution", "surface form": "algebraic"}
-        elif "particle size" in options and options["particle size"] != "distribution":
+        # Necessary/default options
+        options = options or {}
+        if "particle size" in options and options["particle size"] != "distribution":
             raise pybamm.OptionError(
                 "particle size must be 'distribution' for MPM not '{}'".format(
                     options["particle size"]
                 )
             )
-        elif "surface form" in options and options["surface form"] != "algebraic":
+        elif "surface form" in options and options["surface form"] == "false":
             raise pybamm.OptionError(
-                "surface form must be 'algebraic' for MPM not '{}'".format(
-                    options["surface form"]
-                )
+                "surface form must be 'algebraic' or 'differential' for MPM not 'false'"
             )
         else:
-            options["particle size"] = "distribution"
-            options["surface form"] = "algebraic"
+            surface_form = options.get("surface form", "algebraic")
+            options.update(
+                {"particle size": "distribution", "surface form": surface_form}
+            )
         super().__init__(options, name, build)
 
         pybamm.citations.register("Kirk2020")
         pybamm.citations.register("Kirk2021")
 
-    def set_particle_submodel(self):
-        for domain in ["negative", "positive"]:
-            if self.options["particle"] == "Fickian diffusion":
-                submod = pybamm.particle.FickianDiffusion(
-                    self.param, domain, self.options, x_average=True, phase="primary"
-                )
-            elif self.options["particle"] == "uniform profile":
-                submod = pybamm.particle.XAveragedPolynomialProfile(
-                    self.param, domain, self.options, phase="primary"
-                )
-            self.submodels[f"{domain} particle"] = submod
-
     @property
     def default_parameter_values(self):
         default_params = super().default_parameter_values
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
index eb2a88f880..aa8674acca 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -42,6 +42,12 @@ def test_particle_uniform(self):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_differential_surface_form(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.MPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
     def test_voltage_control(self):
         options = {"operating mode": "voltage"}
         model = pybamm.lithium_ion.MPM(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
index 5c824164ea..4e78727e4e 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -55,11 +55,17 @@ def test_particle_quadratic(self):
         with self.assertRaises(NotImplementedError):
             pybamm.lithium_ion.MPM(options)
 
+    def test_differential_surface_form(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
     def test_necessary_options(self):
         options = {"particle size": "single"}
         with self.assertRaises(pybamm.OptionError):
             pybamm.lithium_ion.MPM(options)
-        options = {"surface form": "none"}
+
+        options = {"surface form": "false"}
         with self.assertRaises(pybamm.OptionError):
             pybamm.lithium_ion.MPM(options)