diff --git a/CHANGELOG.md b/CHANGELOG.md
index ea13a5ae37..33953bf8f8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Features
 
 -   Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844))
+-   Added capacitance effects to lithium-ion models ([#842](https://github.com/pybamm-team/PyBaMM/pull/842))
 -   Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824))
 -   Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Added options to export a solution to matlab or csv ([#811](https://github.com/pybamm-team/PyBaMM/pull/811))
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 035e7cb0df..efda022d15 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -224,10 +224,10 @@ def options(self, extra_options):
 
         # Options that are incompatible with models
         if isinstance(self, pybamm.lithium_ion.BaseModel):
-            if options["surface form"] is not False:
-                raise pybamm.OptionError(
-                    "surface form not implemented for lithium-ion models"
-                )
+            # if options["surface form"] is not False:
+            #     raise pybamm.OptionError(
+            #         "surface form not implemented for lithium-ion models"
+            #     )
             if options["convection"] is True:
                 raise pybamm.OptionError(
                     "convection not implemented for lithium-ion models"
diff --git a/pybamm/models/full_battery_models/lithium_ion/dfn.py b/pybamm/models/full_battery_models/lithium_ion/dfn.py
index 0c4bc7c60d..582629c921 100644
--- a/pybamm/models/full_battery_models/lithium_ion/dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/dfn.py
@@ -84,24 +84,40 @@ def set_particle_submodel(self):
 
     def set_solid_submodel(self):
 
-        self.submodels["negative electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Negative", self.reactions
-        )
-        self.submodels["positive electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Positive", self.reactions
-        )
+        if self.options["surface form"] is False:
+            submod_n = pybamm.electrode.ohm.Full(self.param, "Negative", self.reactions)
+            submod_p = pybamm.electrode.ohm.Full(self.param, "Positive", self.reactions)
+        else:
+            submod_n = pybamm.electrode.ohm.SurfaceForm(self.param, "Negative")
+            submod_p = pybamm.electrode.ohm.SurfaceForm(self.param, "Positive")
+
+        self.submodels["negative electrode"] = submod_n
+        self.submodels["positive electrode"] = submod_p
 
     def set_electrolyte_submodel(self):
 
         electrolyte = pybamm.electrolyte.stefan_maxwell
+        surf_form = electrolyte.conductivity.surface_potential_form
 
-        self.submodels["electrolyte conductivity"] = electrolyte.conductivity.Full(
-            self.param, self.reactions
-        )
         self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
             self.param, self.reactions
         )
 
+        if self.options["surface form"] is False:
+            self.submodels["electrolyte conductivity"] = electrolyte.conductivity.Full(
+                self.param, self.reactions
+            )
+        elif self.options["surface form"] == "differential":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    domain.lower() + " electrolyte conductivity"
+                ] = surf_form.FullDifferential(self.param, domain, self.reactions)
+        elif self.options["surface form"] == "algebraic":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    domain.lower() + " electrolyte conductivity"
+                ] = surf_form.FullAlgebraic(self.param, domain, self.reactions)
+
     @property
     def default_geometry(self):
         dimensionality = self.options["dimensionality"]
diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py
index 97e55a2786..023b64bc76 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spm.py
@@ -32,6 +32,7 @@ class SPM(BaseModel):
     def __init__(self, options=None, name="Single Particle Model", build=True):
         super().__init__(options, name)
 
+        self.set_reactions()
         self.set_external_circuit_submodel()
         self.set_porosity_submodel()
         self.set_tortuosity_submodels()
@@ -57,12 +58,21 @@ def set_convection_submodel(self):
 
     def set_interfacial_submodel(self):
 
-        self.submodels[
-            "negative interface"
-        ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Negative")
-        self.submodels[
-            "positive interface"
-        ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Positive")
+        if self.options["surface form"] is False:
+            self.submodels[
+                "negative interface"
+            ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Negative")
+            self.submodels[
+                "positive interface"
+            ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Positive")
+        else:
+            self.submodels[
+                "negative interface"
+            ] = pybamm.interface.lithium_ion.ButlerVolmer(self.param, "Negative")
+
+            self.submodels[
+                "positive interface"
+            ] = pybamm.interface.lithium_ion.ButlerVolmer(self.param, "Positive")
 
     def set_particle_submodel(self):
 
@@ -96,10 +106,26 @@ def set_positive_electrode_submodel(self):
     def set_electrolyte_submodel(self):
 
         electrolyte = pybamm.electrolyte.stefan_maxwell
+        surf_form = electrolyte.conductivity.surface_potential_form
 
-        self.submodels[
-            "electrolyte conductivity"
-        ] = electrolyte.conductivity.LeadingOrder(self.param)
+        if self.options["surface form"] is False:
+            self.submodels[
+                "leading-order electrolyte conductivity"
+            ] = electrolyte.conductivity.LeadingOrder(self.param)
+
+        elif self.options["surface form"] == "differential":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    "leading-order " + domain.lower() + " electrolyte conductivity"
+                ] = surf_form.LeadingOrderDifferential(
+                    self.param, domain, self.reactions
+                )
+
+        elif self.options["surface form"] == "algebraic":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    "leading-order " + domain.lower() + " electrolyte conductivity"
+                ] = surf_form.LeadingOrderAlgebraic(self.param, domain, self.reactions)
         self.submodels[
             "electrolyte diffusion"
         ] = electrolyte.diffusion.ConstantConcentration(self.param)
diff --git a/pybamm/models/submodels/interface/lithium_ion.py b/pybamm/models/submodels/interface/lithium_ion.py
index e995a9eb75..f4cc6b9d88 100644
--- a/pybamm/models/submodels/interface/lithium_ion.py
+++ b/pybamm/models/submodels/interface/lithium_ion.py
@@ -3,6 +3,7 @@
 #
 from .base_interface import BaseInterface
 from . import inverse_kinetics, kinetics
+import pybamm
 
 
 class BaseInterfaceLithiumIon(BaseInterface):
@@ -43,6 +44,16 @@ def _get_exchange_current_density(self, variables):
         c_e = variables[self.domain + " electrolyte concentration"]
         T = variables[self.domain + " electrode temperature"]
 
+        # If variable was broadcast, take only the orphan
+        if (
+            isinstance(c_s_surf, pybamm.Broadcast)
+            and isinstance(c_e, pybamm.Broadcast)
+            and isinstance(T, pybamm.Broadcast)
+        ):
+            c_s_surf = c_s_surf.orphans[0]
+            c_e = c_e.orphans[0]
+            T = T.orphans[0]
+
         if self.domain == "Negative":
             prefactor = self.param.m_n(T) / self.param.C_r_n
 
@@ -75,6 +86,11 @@ def _get_open_circuit_potential(self, variables):
         c_s_surf = variables[self.domain + " particle surface concentration"]
         T = variables[self.domain + " electrode temperature"]
 
+        # If variable was broadcast, take only the orphan
+        if isinstance(c_s_surf, pybamm.Broadcast) and isinstance(T, pybamm.Broadcast):
+            c_s_surf = c_s_surf.orphans[0]
+            T = T.orphans[0]
+
         if self.domain == "Negative":
             ocp = self.param.U_n(c_s_surf, T)
             dUdT = self.param.dUdT_n(c_s_surf)
diff --git a/pybamm/parameters/print_parameters.py b/pybamm/parameters/print_parameters.py
index 381c3f369d..a78b9f31f7 100644
--- a/pybamm/parameters/print_parameters.py
+++ b/pybamm/parameters/print_parameters.py
@@ -67,7 +67,10 @@ def print_parameters(parameters, parameter_values, output_file=None):
                 proc_symbol = parameter_values.process_symbol(symbol)
                 if not (
                     callable(proc_symbol)
-                    or isinstance(proc_symbol, pybamm.Concatenation)
+                    or any(
+                        isinstance(x, (pybamm.Concatenation, pybamm.Broadcast))
+                        for x in proc_symbol.pre_order()
+                    )
                 ):
                     evaluated_parameters[name].append(proc_symbol.evaluate(t=0))
 
diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py
index 865cb80842..6fdd5232fc 100644
--- a/pybamm/parameters/standard_parameters_lithium_ion.py
+++ b/pybamm/parameters/standard_parameters_lithium_ion.py
@@ -92,7 +92,12 @@
 # Electrochemical reactions
 ne_n = pybamm.Parameter("Negative electrode electrons in reaction")
 ne_p = pybamm.Parameter("Positive electrode electrons in reaction")
-C_dl_dimensional = pybamm.Parameter("Double-layer capacity [F.m-2]")
+C_dl_n_dimensional = pybamm.Parameter(
+    "Negative electrode double-layer capacity [F.m-2]"
+)
+C_dl_p_dimensional = pybamm.Parameter(
+    "Positive electrode double-layer capacity [F.m-2]"
+)
 
 
 # Initial conditions
@@ -285,10 +290,15 @@ def U_p_dimensional(sto, T):
 centre_z_tab_p = pybamm.geometric_parameters.centre_z_tab_p
 
 # Microscale geometry
-var = pybamm.standard_spatial_vars
-epsilon_n = pybamm.FunctionParameter("Negative electrode porosity", var.x_n)
-epsilon_s = pybamm.FunctionParameter("Separator porosity", var.x_s)
-epsilon_p = pybamm.FunctionParameter("Positive electrode porosity", var.x_p)
+epsilon_n = pybamm.FunctionParameter(
+    "Negative electrode porosity", pybamm.standard_spatial_vars.x_n
+)
+epsilon_s = pybamm.FunctionParameter(
+    "Separator porosity", pybamm.standard_spatial_vars.x_s
+)
+epsilon_p = pybamm.FunctionParameter(
+    "Positive electrode porosity", pybamm.standard_spatial_vars.x_p
+)
 epsilon = pybamm.Concatenation(epsilon_n, epsilon_s, epsilon_p)
 epsilon_s_n = pybamm.Parameter("Negative electrode active material volume fraction")
 epsilon_s_p = pybamm.Parameter("Positive electrode active material volume fraction")
@@ -315,7 +325,7 @@ def U_p_dimensional(sto, T):
 
 # Electrolyte Properties
 t_plus = pybamm.Parameter("Cation transference number")
-beta_surf = 0
+beta_surf = pybamm.Scalar(0)
 s = 1 - t_plus
 
 
@@ -329,10 +339,10 @@ def chi(c_e):
 
 # Electrochemical Reactions
 C_dl_n = (
-    C_dl_dimensional * potential_scale / interfacial_current_scale_n / tau_discharge
+    C_dl_n_dimensional * potential_scale / interfacial_current_scale_n / tau_discharge
 )
 C_dl_p = (
-    C_dl_dimensional * potential_scale / interfacial_current_scale_p / tau_discharge
+    C_dl_p_dimensional * potential_scale / interfacial_current_scale_p / tau_discharge
 )
 
 # Electrical
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index 96f97f18eb..23d74a461c 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -149,22 +149,24 @@ def initialise_2D(self):
         self.entries = entries
         self.dimensions = 2
         if self.domain[0] in ["negative particle", "positive particle"]:
-            self.spatial_var_name = "r"
+            self.first_dimension = "r"
             self.r_sol = space
         elif self.domain[0] in [
             "negative electrode",
             "separator",
             "positive electrode",
         ]:
-            self.spatial_var_name = "x"
+            self.first_dimension = "x"
             self.x_sol = space
         elif self.domain == ["current collector"]:
-            self.spatial_var_name = "z"
+            self.first_dimension = "z"
             self.z_sol = space
         else:
-            self.spatial_var_name = "x"
+            self.first_dimension = "x"
             self.x_sol = space
 
+        self.first_dim_pts = space
+
         # set up interpolation
         # note that the order of 't' and 'space' is the reverse of what you'd expect
 
@@ -242,6 +244,8 @@ def initialise_3D(self):
         # assign attributes for reference
         self.entries = entries
         self.dimensions = 3
+        self.first_dim_pts = first_dim_pts
+        self.second_dim_pts = second_dim_pts
 
         # set up interpolation
         self._interpolation_function = interp.RegularGridInterpolator(
@@ -335,7 +339,7 @@ def __call__(self, t=None, x=None, r=None, y=None, z=None, warn=True):
 
     def call_2D(self, t, x, r, z):
         "Evaluate a 2D variable"
-        spatial_var = eval_dimension_name(self.spatial_var_name, x, r, None, z)
+        spatial_var = eval_dimension_name(self.first_dimension, x, r, None, z)
         return self._interpolation_function(t, spatial_var)
 
     def call_3D(self, t, x, r, y, z):
diff --git a/pybamm/quick_plot.py b/pybamm/quick_plot.py
index 4d44fb8d77..9c49794e20 100644
--- a/pybamm/quick_plot.py
+++ b/pybamm/quick_plot.py
@@ -195,8 +195,8 @@ def set_output_variables(self, output_variables, solutions):
 
             # Set the x variable for any two-dimensional variables
             if first_variable.dimensions == 2:
-                spatial_variable_key = first_variable.spatial_var_name
-                spatial_variable_value = first_variable.mesh[0].edges
+                spatial_variable_key = first_variable.first_dimension
+                spatial_variable_value = first_variable.first_dim_pts
                 self.spatial_variable[key] = (
                     spatial_variable_key,
                     spatial_variable_value,
diff --git a/tests/integration/test_models/standard_output_comparison.py b/tests/integration/test_models/standard_output_comparison.py
index e4aa69e0f6..b36ac1a7bc 100644
--- a/tests/integration/test_models/standard_output_comparison.py
+++ b/tests/integration/test_models/standard_output_comparison.py
@@ -51,9 +51,9 @@ def test_all(self, skip_first_timestep=False):
         self.run_test_class(VariablesComparison, skip_first_timestep)
 
         if self.chemistry == "Lithium-ion":
-            self.run_test_class(ParticleConcentrationComparison)
+            self.run_test_class(ParticleConcentrationComparison, skip_first_timestep)
         elif self.chemistry == "Lead-acid":
-            self.run_test_class(PorosityComparison)
+            self.run_test_class(PorosityComparison, skip_first_timestep)
 
 
 class BaseOutputComparison(object):
@@ -67,13 +67,14 @@ def compare(self, var, tol=1e-2):
         model_variables = [solution[var] for solution in self.solutions]
         var0 = model_variables[0]
 
-        if var0.mesh is None:
-            x = None
-        else:
-            x = var0.mesh[0].nodes
+        spatial_pts = {}
+        if var0.dimensions >= 2:
+            spatial_pts[var0.first_dimension] = var0.first_dim_pts
+        if var0.dimensions >= 3:
+            spatial_pts[var0.second_dimension] = var0.second_dim_pts
 
         # Calculate tolerance based on the value of var0
-        maxvar0 = np.max(abs(var0(self.t, x)))
+        maxvar0 = np.max(abs(var0(self.t, **spatial_pts)))
         if maxvar0 < 1e-14:
             decimal = -int(np.log10(tol))
         else:
@@ -82,7 +83,7 @@ def compare(self, var, tol=1e-2):
         for model_var in model_variables[1:]:
             np.testing.assert_equal(var0.dimensions, model_var.dimensions)
             np.testing.assert_array_almost_equal(
-                model_var(self.t, x), var0(self.t, x), decimal
+                model_var(self.t, **spatial_pts), var0(self.t, **spatial_pts), decimal
             )
 
 
@@ -123,7 +124,6 @@ def test_all(self):
         self.compare("X-averaged negative electrode open circuit potential")
         self.compare("X-averaged positive electrode open circuit potential")
         self.compare("Terminal voltage")
-        self.compare("X-averaged electrolyte overpotential")
         self.compare("X-averaged solid phase ohmic losses")
         self.compare("Negative electrode reaction overpotential")
         self.compare("Positive electrode reaction overpotential")
diff --git a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py
index a3700d62d1..30b2e59804 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py
@@ -48,10 +48,9 @@ def test_compare_averages_asymptotics(self):
         comparison = StandardOutputComparison(solutions)
         comparison.test_averages()
 
-    def test_compare_outputs_capacitance(self):
+    def test_compare_outputs_surface_form(self):
         """
-        Check that the leading-order model solution converges linearly in C_e to the
-        full model solution
+        Check that the models agree with the different surface forms
         """
         # load models
         options = [
@@ -59,6 +58,7 @@ def test_compare_outputs_capacitance(self):
         ]
         model_combos = [
             ([pybamm.lead_acid.LOQS(opt) for opt in options]),
+            ([pybamm.lead_acid.Composite(opt) for opt in options]),
             ([pybamm.lead_acid.Full(opt) for opt in options]),
         ]
 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs.py
new file mode 100644
index 0000000000..973579d628
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs.py
@@ -0,0 +1,61 @@
+#
+# Tests for the surface formulation
+#
+import pybamm
+import numpy as np
+import unittest
+from tests import StandardOutputComparison
+
+
+class TestCompareOutputs(unittest.TestCase):
+    def test_compare_outputs_surface_form(self):
+        # load models
+        options = [
+            {"surface form": cap} for cap in [False, "differential", "algebraic"]
+        ]
+        model_combos = [
+            ([pybamm.lithium_ion.SPM(opt) for opt in options]),
+            ([pybamm.lithium_ion.SPMe(opt) for opt in options]),
+            ([pybamm.lithium_ion.DFN(opt) for opt in options]),
+        ]
+
+        for models in model_combos:
+            # load parameter values (same for all models)
+            param = models[0].default_parameter_values
+            param.update({"Current function [A]": 1})
+            for model in models:
+                param.process_model(model)
+
+            # set mesh
+            var = pybamm.standard_spatial_vars
+            var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5}
+
+            # discretise models
+            discs = {}
+            for model in models:
+                geometry = model.default_geometry
+                param.process_geometry(geometry)
+                mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+                disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+                disc.process_model(model)
+                discs[model] = disc
+
+            # solve model
+            solutions = []
+            t_eval = np.linspace(0, 3600, 100)
+            for model in models:
+                solution = pybamm.CasadiSolver().solve(model, t_eval)
+                solutions.append(solution)
+
+            # compare outputs
+            comparison = StandardOutputComparison(solutions)
+            comparison.test_all(skip_first_timestep=True)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 8082aae7d2..be6eb491f8 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -95,6 +95,18 @@ def test_particle_fast_diffusion(self):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
     def test_particle_distribution_in_x(self):
         model = pybamm.lithium_ion.DFN()
         param = model.default_parameter_values
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index cbed675bf2..9a933dade6 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -103,6 +103,18 @@ def test_particle_fast_diffusion(self):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index bf48669845..0709b96686 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -95,6 +95,18 @@ def test_particle_fast_diffusion(self):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPMe(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPMe(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index 901f43f1c0..511b4c3541 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -13,7 +13,7 @@ def test_set_up(self):
                 "Charge at 1 A until 4.1 V",
                 "Hold at 4.1 V until 50 mA",
                 "Discharge at 2 W for 1 hour",
-            ],
+            ]
         )
         model = pybamm.lithium_ion.DFN()
         sim = pybamm.Simulation(model, experiment=experiment)
@@ -48,7 +48,7 @@ def test_set_up(self):
         self.assertEqual(sim._experiment_inputs[3]["Voltage cut-off [V]"], -1e10)
 
         self.assertEqual(
-            sim._experiment_times, [3600, 7 * 24 * 3600, 7 * 24 * 3600, 3600],
+            sim._experiment_times, [3600, 7 * 24 * 3600, 7 * 24 * 3600, 3600]
         )
 
         self.assertIn(
@@ -75,11 +75,11 @@ def test_run_experiment(self):
                 "Charge at 1 A until 4.1 V",
                 "Hold at 4.1 V until C/2",
                 "Discharge at 2 W for 1 hour",
-            ],
+            ]
         )
         model = pybamm.lithium_ion.SPM()
         sim = pybamm.Simulation(model, experiment=experiment)
-        sim.solve()
+        sim.solve(solver=pybamm.CasadiSolver())
         self.assertEqual(sim._solution.termination, "final time")
 
     def test_run_experiment_breaks_early(self):
@@ -89,7 +89,7 @@ def test_run_experiment_breaks_early(self):
         pybamm.set_logging_level("ERROR")
         # giving the time, should get ignored
         t_eval = [0, 1]
-        sim.solve(t_eval)
+        sim.solve(t_eval, solver=pybamm.CasadiSolver())
         pybamm.set_logging_level("WARNING")
         self.assertIn("event", sim._solution.termination)
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py
index 183fe3cc46..500a00d37f 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py
@@ -54,11 +54,6 @@ def test_well_posed(self):
 
 
 class TestLeadAcidCompositeWithSideReactions(unittest.TestCase):
-    def test_well_posed_differential(self):
-        options = {"surface form": "differential", "side reactions": ["oxygen"]}
-        model = pybamm.lead_acid.Composite(options)
-        model.check_well_posedness()
-
     def test_well_posed_algebraic(self):
         options = {"surface form": "algebraic", "side reactions": ["oxygen"]}
         model = pybamm.lead_acid.Composite(options)
@@ -69,10 +64,6 @@ def test_well_posed_algebraic(self):
 
 
 class TestLeadAcidCompositeExtended(unittest.TestCase):
-    def test_well_posed(self):
-        model = pybamm.lead_acid.CompositeExtended()
-        model.check_well_posedness()
-
     def test_well_posed_differential_side_reactions(self):
         options = {"surface form": "differential", "side reactions": ["oxygen"]}
         model = pybamm.lead_acid.CompositeExtended(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
index 3ac309316f..9e5bda3d5b 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
@@ -7,8 +7,6 @@
 
 class TestBaseLithiumIonModel(unittest.TestCase):
     def test_incompatible_options(self):
-        with self.assertRaisesRegex(pybamm.OptionError, "surface form not implemented"):
-            pybamm.lithium_ion.BaseModel({"surface form": "differential"})
         with self.assertRaisesRegex(pybamm.OptionError, "convection not implemented"):
             pybamm.lithium_ion.BaseModel({"convection": True})
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 4fcabbccb0..f6517319dc 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -164,6 +164,16 @@ def test_particle_fast_diffusion(self):
         model = pybamm.lithium_ion.DFN(options)
         model.check_well_posedness()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.DFN(options)
+        model.check_well_posedness()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.DFN(options)
+        model.check_well_posedness()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 87da7b98b2..9d75316b73 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -177,6 +177,16 @@ def test_particle_fast_diffusion(self):
         model = pybamm.lithium_ion.SPM(options)
         model.check_well_posedness()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPM(options)
+        model.check_well_posedness()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPM(options)
+        model.check_well_posedness()
+
 
 class TestSPMExternalCircuits(unittest.TestCase):
     def test_well_posed_voltage(self):
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index fb1e48dac7..9020d70b53 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -164,6 +164,16 @@ def test_particle_fast_diffusion(self):
         model = pybamm.lithium_ion.SPMe(options)
         model.check_well_posedness()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPMe(options)
+        model.check_well_posedness()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPMe(options)
+        model.check_well_posedness()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")