diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3e743a4b5..8af4eba2a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#18](https://github.com/pybop-team/PyBOP/pull/18) - Adds geometric parameter fitting capability, via `model.rebuild()` with `model.rebuild_parameters`.
 - [#203](https://github.com/pybop-team/PyBOP/pull/203) - Adds support for modern Python packaging via a `pyproject.toml` file and configures the `pytest` test runner and `ruff` linter to use their configurations stored as declarative metadata.
 - [#123](https://github.com/pybop-team/PyBOP/issues/123) - Configures scheduled tests to run against the last three PyPI releases of PyBaMM via dynamic GitHub Actions matrix generation.
 - [#187](https://github.com/pybop-team/PyBOP/issues/187) - Adds M1 Github runner to `test_on_push` workflow, updt. self-hosted supported python versions in scheduled tests.
diff --git a/examples/scripts/spm_CMAES.py b/examples/scripts/spm_CMAES.py
index e24eb7a15..f5728019e 100644
--- a/examples/scripts/spm_CMAES.py
+++ b/examples/scripts/spm_CMAES.py
@@ -8,14 +8,16 @@
 # Fitting parameters
 parameters = [
     pybop.Parameter(
-        "Negative electrode active material volume fraction",
-        prior=pybop.Gaussian(0.6, 0.05),
-        bounds=[0.5, 0.8],
+        "Negative particle radius [m]",
+        prior=pybop.Gaussian(6e-06, 0.1e-6),
+        bounds=[1e-6, 9e-6],
+        true_value=parameter_set["Negative particle radius [m]"],
     ),
     pybop.Parameter(
-        "Positive electrode active material volume fraction",
-        prior=pybop.Gaussian(0.48, 0.05),
-        bounds=[0.4, 0.7],
+        "Positive particle radius [m]",
+        prior=pybop.Gaussian(4.5e-06, 0.1e-6),
+        bounds=[1e-6, 9e-6],
+        true_value=parameter_set["Positive particle radius [m]"],
     ),
 ]
 
@@ -42,6 +44,13 @@
 
 # Run the optimisation
 x, final_cost = optim.run()
+print(
+    "True parameters:",
+    [
+        parameters[0].true_value,
+        parameters[1].true_value,
+    ],
+)
 print("Estimated parameters:", x)
 
 # Plot the timeseries output
@@ -57,5 +66,4 @@
 pybop.plot_cost2d(cost, steps=15)
 
 # Plot the cost landscape with optimisation path and updated bounds
-bounds = np.array([[0.55, 0.75], [0.48, 0.68]])
-pybop.plot_cost2d(cost, optim=optim, bounds=bounds, steps=15)
+pybop.plot_cost2d(cost, optim=optim, steps=15)
diff --git a/examples/standalone/model.py b/examples/standalone/model.py
index 9943ff3d8..2295d080a 100644
--- a/examples/standalone/model.py
+++ b/examples/standalone/model.py
@@ -56,3 +56,4 @@ def __init__(
         self._built_initial_soc = None
         self._mesh = None
         self._disc = None
+        self.rebuild_parameters = {}
diff --git a/pybop/_problem.py b/pybop/_problem.py
index 87b933b04..38e28ce20 100644
--- a/pybop/_problem.py
+++ b/pybop/_problem.py
@@ -61,7 +61,7 @@ def __init__(
 
         # Add the initial values to the parameter definitions
         for i, param in enumerate(self.parameters):
-            param.update(value=self.x0[i])
+            param.update(initial_value=self.x0[i])
 
     def evaluate(self, x):
         """
@@ -118,6 +118,10 @@ def target(self):
         """
         return self._target
 
+    @property
+    def model(self):
+        return self._model
+
 
 class FittingProblem(BaseProblem):
     """
@@ -149,6 +153,7 @@ def __init__(
     ):
         super().__init__(parameters, model, check_model, signal, init_soc, x0)
         self._dataset = dataset.data
+        self.x = self.x0
 
         # Check that the dataset contains time and current
         for name in ["Time [s]", "Current function [A]"] + self.signal:
@@ -199,6 +204,13 @@ def evaluate(self, x):
         y : np.ndarray
             The model output y(t) simulated with inputs x.
         """
+        if (x != self.x).any() and self._model.matched_parameters:
+            for i, param in enumerate(self.parameters):
+                param.update(value=x[i])
+
+            self._model.rebuild(parameters=self.parameters)
+            self.x = x
+
         y = np.asarray(self._model.simulate(inputs=x, t_eval=self._time_data))
 
         return y
@@ -218,6 +230,10 @@ def evaluateS1(self, x):
             A tuple containing the simulation result y(t) and the sensitivities dy/dx(t) evaluated
             with given inputs x.
         """
+        if self._model.matched_parameters:
+            raise RuntimeError(
+                "Gradient not available when using geometric parameters."
+            )
 
         y, dy = self._model.simulateS1(
             inputs=x,
diff --git a/pybop/costs/design_costs.py b/pybop/costs/design_costs.py
index 297697e39..8b02f1f28 100644
--- a/pybop/costs/design_costs.py
+++ b/pybop/costs/design_costs.py
@@ -41,7 +41,7 @@ def __init__(self, problem, update_capacity=False):
             )
         warnings.warn(nominal_capacity_warning, UserWarning)
         self.update_capacity = update_capacity
-        self.parameter_set = problem._model._parameter_set
+        self.parameter_set = problem.model.parameter_set
         self.update_simulation_data(problem.x0)
 
     def update_simulation_data(self, initial_conditions):
@@ -54,7 +54,7 @@ def update_simulation_data(self, initial_conditions):
             The initial conditions for the simulation.
         """
         if self.update_capacity:
-            self.problem._model.approximate_capacity(self.problem.x0)
+            self.problem.model.approximate_capacity(self.problem.x0)
         solution = self.problem.evaluate(initial_conditions)
         self.problem._time_data = solution[:, -1]
         self.problem._target = solution[:, 0:-1]
@@ -116,12 +116,12 @@ def _evaluate(self, x, grad=None):
                 warnings.filterwarnings("error", category=UserWarning)
 
                 if self.update_capacity:
-                    self.problem._model.approximate_capacity(x)
+                    self.problem.model.approximate_capacity(x)
                 solution = self.problem.evaluate(x)
 
                 voltage, current = solution[:, 0], solution[:, 1]
                 negative_energy_density = -np.trapz(voltage * current, dx=self.dt) / (
-                    3600 * self.problem._model.cell_mass(self.parameter_set)
+                    3600 * self.problem.model.cell_mass(self.parameter_set)
                 )
 
                 return negative_energy_density
@@ -166,12 +166,12 @@ def _evaluate(self, x, grad=None):
                 warnings.filterwarnings("error", category=UserWarning)
 
                 if self.update_capacity:
-                    self.problem._model.approximate_capacity(x)
+                    self.problem.model.approximate_capacity(x)
                 solution = self.problem.evaluate(x)
 
                 voltage, current = solution[:, 0], solution[:, 1]
                 negative_energy_density = -np.trapz(voltage * current, dx=self.dt) / (
-                    3600 * self.problem._model.cell_volume(self.parameter_set)
+                    3600 * self.problem.model.cell_volume(self.parameter_set)
                 )
 
                 return negative_energy_density
diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py
index a02aeca6e..2da337f18 100644
--- a/pybop/costs/fitting_costs.py
+++ b/pybop/costs/fitting_costs.py
@@ -37,7 +37,6 @@ def _evaluate(self, x, grad=None):
             The root mean square error.
 
         """
-
         prediction = self.problem.evaluate(x)
 
         if len(prediction) < len(self._target):
diff --git a/pybop/models/base_model.py b/pybop/models/base_model.py
index 945ee1bdb..1f56f9b4e 100644
--- a/pybop/models/base_model.py
+++ b/pybop/models/base_model.py
@@ -2,6 +2,7 @@
 from dataclasses import dataclass
 from typing import Any, Dict, Optional
 import pybamm
+import copy
 import numpy as np
 import casadi
 
@@ -57,6 +58,9 @@ def __init__(self, name="Base Model"):
         self.parameters = None
         self.dataset = None
         self.signal = None
+        self.matched_parameters = {}
+        self.non_matched_parameters = {}
+        self.fit_keys = []
         self.param_check_counter = 0
         self.allow_infeasible_solutions = True
 
@@ -87,10 +91,11 @@ def build(
         """
         self.dataset = dataset
         self.parameters = parameters
-        if self.parameters is None:
-            self.fit_keys = []
-        else:
+        if self.parameters is not None:
+            self.set_parameter_classification(self.parameters)
             self.fit_keys = [param.name for param in self.parameters]
+        else:
+            self.fit_keys = []
 
         if init_soc is not None:
             self.set_init_soc(init_soc)
@@ -140,24 +145,24 @@ def set_init_soc(self, init_soc):
         # Save solved initial SOC in case we need to rebuild the model
         self._built_initial_soc = init_soc
 
-    def set_params(self):
+    def set_params(self, rebuild=False):
         """
         Assign the parameters to the model.
 
         This method processes the model with the given parameters, sets up
         the geometry, and updates the model instance.
         """
-        if self.model_with_set_params:
+        if self.model_with_set_params and not rebuild:
             return
 
         # Mark any simulation inputs in the parameter set
-        if self.parameters is not None:
+        if self.non_matched_parameters:
             for i in self.fit_keys:
                 self._parameter_set[i] = "[input]"
 
-        if self.dataset is not None and self.parameters is not None:
+        if self.dataset is not None and self.non_matched_parameters:
             if "Current function [A]" not in self.fit_keys:
-                self.parameter_set["Current function [A]"] = pybamm.Interpolant(
+                self._parameter_set["Current function [A]"] = pybamm.Interpolant(
                     self.dataset["Time [s]"],
                     self.dataset["Current function [A]"],
                     pybamm.t,
@@ -172,6 +177,93 @@ def set_params(self):
             self._parameter_set.process_geometry(self.geometry)
         self.pybamm_model = self._model_with_set_params
 
+    def rebuild(
+        self,
+        dataset=None,
+        parameters=None,
+        parameter_set=None,
+        check_model=True,
+        init_soc=None,
+    ):
+        """
+        Rebuild the PyBaMM model for a given parameter set.
+
+        This method requires the self.build() method to be called first, and
+        then rebuilds the model for a given parameter set. Specifically,
+        this method applies the given parameters, sets up the mesh and discretization if needed, and prepares the model
+        for simulations.
+
+        Parameters
+        ----------
+        dataset : pybamm.Dataset, optional
+            The dataset to be used in the model construction.
+        parameters : dict, optional
+            A dictionary containing parameter values to apply to the model.
+        parameter_set : pybop.parameter_set, optional
+            A PyBOP parameter set object or a dictionary containing the parameter values
+        check_model : bool, optional
+            If True, the model will be checked for correctness after construction.
+        init_soc : float, optional
+            The initial state of charge to be used in simulations.
+        """
+        self.dataset = dataset
+        self.parameters = parameters
+        if parameters is not None:
+            self.set_parameter_classification(parameters)
+
+        if init_soc is not None:
+            self.set_init_soc(init_soc)
+
+        if self._built_model is None:
+            raise ValueError("Model must be built before calling rebuild")
+
+        self.set_params(rebuild=True)
+        self._mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.var_pts)
+        self._disc = pybamm.Discretisation(self.mesh, self.spatial_methods)
+        self._built_model = self._disc.process_model(
+            self._model_with_set_params, inplace=False, check_model=check_model
+        )
+
+        # Clear solver and setup model
+        self._solver._model_set_up = {}
+
+    def set_parameter_classification(self, parameters):
+        """
+        Set the parameter classification for the model.
+
+        Parameters
+        ----------
+        parameters : Pybop.ParameterSet
+
+        Returns
+        -------
+        None
+            The method updates attributes on self.
+
+        """
+        processed_parameters = {param.name: param.value for param in parameters}
+        matched_parameters = {
+            param: processed_parameters[param]
+            for param in processed_parameters
+            if param in self.rebuild_parameters
+        }
+        non_matched_parameters = {
+            param: processed_parameters[param]
+            for param in processed_parameters
+            if param not in self.rebuild_parameters
+        }
+
+        self.matched_parameters.update(matched_parameters)
+        self.non_matched_parameters.update(non_matched_parameters)
+
+        if self.matched_parameters:
+            self._parameter_set.update(self.matched_parameters)
+            self._unprocessed_parameter_set = self._parameter_set
+            self.geometry = self.pybamm_model.default_geometry
+
+        if self.non_matched_parameters:
+            self.fit_keys = list(self.non_matched_parameters.keys())
+
     def reinit(
         self, inputs: Inputs, t: float = 0.0, x: Optional[np.ndarray] = None
     ) -> TimeSeriesState:
@@ -243,25 +335,29 @@ def simulate(self, inputs, t_eval) -> np.ndarray[np.float64]:
         ValueError
             If the model has not been built before simulation.
         """
-
         if self._built_model is None:
             raise ValueError("Model must be built before calling simulate")
         else:
-            if not isinstance(inputs, dict):
-                inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)}
+            if not self.fit_keys and self.matched_parameters:
+                sol = self.solver.solve(self.built_model, t_eval=t_eval)
 
-            if self.check_params(
-                inputs=inputs,
-                allow_infeasible_solutions=self.allow_infeasible_solutions,
-            ):
-                sol = self._solver.solve(self.built_model, inputs=inputs, t_eval=t_eval)
+            else:
+                if not isinstance(inputs, dict):
+                    inputs = {key: inputs[i] for i, key in enumerate(self.fit_keys)}
 
-                predictions = [sol[signal].data for signal in self.signal]
+                if self.check_params(
+                    inputs=inputs,
+                    allow_infeasible_solutions=self.allow_infeasible_solutions,
+                ):
+                    sol = self.solver.solve(
+                        self.built_model, inputs=inputs, t_eval=t_eval
+                    )
+                else:
+                    return [np.inf]
 
-                return np.vstack(predictions).T
+            predictions = [sol[signal].data for signal in self.signal]
 
-            else:
-                return [np.inf]
+            return np.vstack(predictions).T
 
     def simulateS1(self, inputs, t_eval):
         """
@@ -449,6 +545,17 @@ def _check_params(self, inputs=None, allow_infeasible_solutions=True):
         """
         return True
 
+    def copy(self):
+        """
+        Return a copy of the model.
+
+        Returns
+        -------
+        BaseModel
+            A copy of the model.
+        """
+        return copy.copy(self)
+
     def cell_mass(self, parameter_set=None):
         """
         Calculate the cell mass in kilograms.
diff --git a/pybop/models/empirical/ecm.py b/pybop/models/empirical/ecm.py
index 501a10494..6c0427c8c 100644
--- a/pybop/models/empirical/ecm.py
+++ b/pybop/models/empirical/ecm.py
@@ -75,6 +75,7 @@ def __init__(
         self._built_initial_soc = None
         self._mesh = None
         self._disc = None
+        self.rebuild_parameters = {}
 
     def _check_params(self, inputs=None, allow_infeasible_solutions=True):
         """
diff --git a/pybop/models/lithium_ion/echem.py b/pybop/models/lithium_ion/echem.py
index afc4a3b2a..bb03f6585 100644
--- a/pybop/models/lithium_ion/echem.py
+++ b/pybop/models/lithium_ion/echem.py
@@ -69,6 +69,7 @@ def __init__(
         self._disc = None
 
         self._electrode_soh = pybamm.lithium_ion.electrode_soh
+        self.rebuild_parameters = self.set_rebuild_parameters()
 
 
 class SPMe(EChemBaseModel):
@@ -140,3 +141,4 @@ def __init__(
         self._disc = None
 
         self._electrode_soh = pybamm.lithium_ion.electrode_soh
+        self.rebuild_parameters = self.set_rebuild_parameters()
diff --git a/pybop/models/lithium_ion/echem_base.py b/pybop/models/lithium_ion/echem_base.py
index 7e9aaabcb..7e5c869fa 100644
--- a/pybop/models/lithium_ion/echem_base.py
+++ b/pybop/models/lithium_ion/echem_base.py
@@ -232,3 +232,28 @@ def approximate_capacity(self, x):
         self._parameter_set.update(
             {"Nominal cell capacity [A.h]": theoretical_capacity}
         )
+
+    def set_rebuild_parameters(self):
+        """
+        Sets the parameters that can be changed when rebuilding the model.
+
+        Returns
+        -------
+        dict
+            A dictionary of parameters that can be changed when rebuilding the model.
+
+        """
+        rebuild_parameters = dict.fromkeys(
+            [
+                "Negative particle radius [m]",
+                "Negative electrode porosity",
+                "Negative electrode thickness [m]",
+                "Positive particle radius [m]",
+                "Positive electrode porosity",
+                "Positive electrode thickness [m]",
+                "Separator porosity",
+                "Separator thickness [m]",
+            ]
+        )
+
+        return rebuild_parameters
diff --git a/pybop/parameters/parameter.py b/pybop/parameters/parameter.py
index 15dffecb7..6136b0acb 100644
--- a/pybop/parameters/parameter.py
+++ b/pybop/parameters/parameter.py
@@ -28,12 +28,15 @@ class Parameter:
         the margin is set outside the interval (0, 1).
     """
 
-    def __init__(self, name, initial_value=None, prior=None, bounds=None):
+    def __init__(
+        self, name, initial_value=None, true_value=None, prior=None, bounds=None
+    ):
         """
         Construct the parameter class with a name, initial value, prior, and bounds.
         """
         self.name = name
         self.prior = prior
+        self.true_value = true_value
         self.initial_value = initial_value
         self.value = initial_value
         self.bounds = bounds
@@ -69,7 +72,7 @@ def rvs(self, n_samples):
 
         return samples
 
-    def update(self, value):
+    def update(self, value=None, initial_value=None):
         """
         Update the parameter's current value.
 
@@ -78,7 +81,12 @@ def update(self, value):
         value : float
             The new value to be assigned to the parameter.
         """
-        self.value = value
+        if value is not None:
+            self.value = value
+        elif initial_value is not None:
+            self.value = initial_value
+        else:
+            raise ValueError("No value provided to update parameter")
 
     def __repr__(self):
         """
diff --git a/pybop/plotting/plot_cost2d.py b/pybop/plotting/plot_cost2d.py
index ee43d9d09..1565ed243 100644
--- a/pybop/plotting/plot_cost2d.py
+++ b/pybop/plotting/plot_cost2d.py
@@ -46,7 +46,7 @@ def plot_cost2d(cost, bounds=None, optim=None, steps=10):
     # Populate cost matrix
     for i, xi in enumerate(x):
         for j, yj in enumerate(y):
-            costs[j, i] = cost([xi, yj])
+            costs[j, i] = cost(np.array([xi, yj]))
 
     # Create figure
     fig = create_figure(x, y, costs, bounds, cost.problem.parameters, optim)
diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py
index 2303f17af..dfce88e48 100644
--- a/tests/unit/test_models.py
+++ b/tests/unit/test_models.py
@@ -77,6 +77,101 @@ def test_build(self):
         model.build()
         assert model.built_model is not None
 
+    @pytest.mark.unit
+    def test_rebuild(self):
+        model = pybop.lithium_ion.SPM()
+        model.build()
+        initial_built_model = model._built_model
+        assert model._built_model is not None
+
+        # Test that the model can be built again
+        model.rebuild()
+        rebuilt_model = model._built_model
+        assert rebuilt_model is not None
+
+        # Filter out special and private attributes
+        attributes_to_compare = [
+            "algebraic",
+            "bcs",
+            "boundary_conditions",
+            "mass_matrix",
+            "parameters",
+            "submodels",
+            "summary_variables",
+            "rhs",
+            "variables",
+            "y_slices",
+        ]
+
+        # Loop through the filtered attributes and compare them
+        for attribute in attributes_to_compare:
+            assert getattr(rebuilt_model, attribute) == getattr(
+                initial_built_model, attribute
+            )
+
+    @pytest.mark.unit
+    def test_rebuild_geometric_parameters(self):
+        parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+        parameters = [
+            pybop.Parameter(
+                "Positive particle radius [m]",
+                prior=pybop.Gaussian(4.8e-06, 0.05e-06),
+                bounds=[4e-06, 6e-06],
+                initial_value=4.8e-06,
+            ),
+            pybop.Parameter(
+                "Negative electrode thickness [m]",
+                prior=pybop.Gaussian(40e-06, 1e-06),
+                bounds=[30e-06, 50e-06],
+                initial_value=48e-06,
+            ),
+        ]
+
+        model = pybop.lithium_ion.SPM(parameter_set=parameter_set)
+        model.build(parameters=parameters)
+        initial_built_model = model.copy()
+        assert initial_built_model._built_model is not None
+
+        # Run prediction
+        t_eval = np.linspace(0, 100, 100)
+        out_init = initial_built_model.predict(t_eval=t_eval)
+
+        # Test that the model can be rebuilt with different geometric parameters
+        parameters[0].update(5e-06)
+        parameters[1].update(45e-06)
+        model.rebuild(parameters=parameters)
+        rebuilt_model = model
+        assert rebuilt_model._built_model is not None
+
+        # Test model geometry
+        assert (
+            rebuilt_model._mesh["negative electrode"].nodes[1]
+            != initial_built_model._mesh["negative electrode"].nodes[1]
+        )
+        assert (
+            rebuilt_model.geometry["negative electrode"]["x_n"]["max"]
+            != initial_built_model.geometry["negative electrode"]["x_n"]["max"]
+        )
+
+        assert (
+            rebuilt_model.geometry["positive particle"]["r_p"]["max"]
+            != initial_built_model.geometry["positive particle"]["r_p"]["max"]
+        )
+
+        assert (
+            rebuilt_model._mesh["positive particle"].nodes[1]
+            != initial_built_model._mesh["positive particle"].nodes[1]
+        )
+
+        # Compare model results
+        out_rebuild = rebuilt_model.predict(t_eval=t_eval)
+        with pytest.raises(AssertionError):
+            np.testing.assert_allclose(
+                out_init["Terminal voltage [V]"].data,
+                out_rebuild["Terminal voltage [V]"].data,
+                atol=1e-5,
+            )
+
     @pytest.mark.unit
     def test_reinit(self):
         k = 0.1
diff --git a/tests/unit/test_parameters.py b/tests/unit/test_parameters.py
new file mode 100644
index 000000000..52f0c69fd
--- /dev/null
+++ b/tests/unit/test_parameters.py
@@ -0,0 +1,53 @@
+import pybop
+import pytest
+
+
+class TestParameters:
+    """
+    A class to test the parameter classes.
+    """
+
+    @pytest.fixture
+    def parameter(self):
+        return pybop.Parameter(
+            "Negative electrode active material volume fraction",
+            prior=pybop.Gaussian(0.6, 0.02),
+            bounds=[0.375, 0.7],
+            initial_value=0.6,
+        )
+
+    @pytest.mark.unit
+    def test_parameter_construction(self, parameter):
+        assert parameter.name == "Negative electrode active material volume fraction"
+        assert parameter.prior.mean == 0.6
+        assert parameter.prior.sigma == 0.02
+        assert parameter.bounds == [0.375, 0.7]
+        assert parameter.initial_value == 0.6
+
+    @pytest.mark.unit
+    def test_parameter_repr(self, parameter):
+        assert (
+            repr(parameter)
+            == "Parameter: Negative electrode active material volume fraction \n Prior: Gaussian, mean: 0.6, sigma: 0.02 \n Bounds: [0.375, 0.7] \n Value: 0.6"
+        )
+
+    @pytest.mark.unit
+    def test_parameter_rvs(self, parameter):
+        samples = parameter.rvs(n_samples=500)
+        assert (samples >= 0.375).all() and (samples <= 0.7).all()
+
+    @pytest.mark.unit
+    def test_parameter_update(self, parameter):
+        # Test value update
+        parameter.update(value=0.534)
+        assert parameter.value == 0.534
+
+        # Test initial value update
+        parameter.update(initial_value=0.654)
+        assert parameter.value == 0.654
+
+    @pytest.mark.unit
+    def test_parameter_margin(self, parameter):
+        assert parameter.margin == 1e-4
+        parameter.set_margin(margin=1e-3)
+        assert parameter.margin == 1e-3
diff --git a/tests/unit/test_problem.py b/tests/unit/test_problem.py
index b92b89a3a..0cdfd8bce 100644
--- a/tests/unit/test_problem.py
+++ b/tests/unit/test_problem.py
@@ -16,14 +16,14 @@ def model(self):
     def parameters(self):
         return [
             pybop.Parameter(
-                "Negative electrode active material volume fraction",
-                prior=pybop.Gaussian(0.5, 0.02),
-                bounds=[0.375, 0.625],
+                "Negative particle radius [m]",
+                prior=pybop.Gaussian(2e-05, 0.1e-5),
+                bounds=[1e-6, 5e-5],
             ),
             pybop.Parameter(
-                "Positive electrode active material volume fraction",
-                prior=pybop.Gaussian(0.65, 0.02),
-                bounds=[0.525, 0.75],
+                "Positive particle radius [m]",
+                prior=pybop.Gaussian(0.5e-05, 0.1e-5),
+                bounds=[1e-6, 5e-5],
             ),
         ]
 
@@ -44,11 +44,11 @@ def experiment(self):
     @pytest.fixture
     def dataset(self, model, experiment):
         model.parameter_set = model.pybamm_model.default_parameter_values
-        x0 = np.array([0.52, 0.63])
+        x0 = np.array([2e-5, 0.5e-5])
         model.parameter_set.update(
             {
-                "Negative electrode active material volume fraction": x0[0],
-                "Positive electrode active material volume fraction": x0[1],
+                "Negative particle radius [m]": x0[0],
+                "Positive particle radius [m]": x0[1],
             }
         )
         solution = model.predict(experiment=experiment)
@@ -76,12 +76,12 @@ def test_base_problem(self, parameters, model):
         assert problem._model == model
 
         with pytest.raises(NotImplementedError):
-            problem.evaluate([0.5, 0.5])
+            problem.evaluate([1e-5, 1e-5])
         with pytest.raises(NotImplementedError):
-            problem.evaluateS1([0.5, 0.5])
+            problem.evaluateS1([1e-5, 1e-5])
 
         with pytest.raises(ValueError):
-            pybop._problem.BaseProblem(parameters, model=model, signal=[0.5, 0.5])
+            pybop._problem.BaseProblem(parameters, model=model, signal=[1e-5, 1e-5])
 
     @pytest.mark.unit
     def test_fitting_problem(self, parameters, dataset, model, signal):
@@ -98,7 +98,7 @@ def test_fitting_problem(self, parameters, dataset, model, signal):
         assert problem._model._built_model is not None
 
         # Test model.simulate
-        model.simulate(inputs=[0.5, 0.5], t_eval=np.linspace(0, 10, 100))
+        model.simulate(inputs=[1e-5, 1e-5], t_eval=np.linspace(0, 10, 100))
 
         # Test problem construction errors
         for bad_dataset in [
@@ -147,5 +147,27 @@ def test_design_problem(self, parameters, experiment, model):
         )  # building postponed with input experiment
 
         # Test model.predict
-        model.predict(inputs=[0.5, 0.5], experiment=experiment)
-        model.predict(inputs=[1.1, 0.5], experiment=experiment)
+        model.predict(inputs=[1e-5, 1e-5], experiment=experiment)
+        model.predict(inputs=[3e-5, 3e-5], experiment=experiment)
+
+    @pytest.mark.unit
+    def test_problem_construct_with_model_predict(
+        self, parameters, model, dataset, signal
+    ):
+        # Construct model and predict
+        out = model.predict(inputs=[1e-5, 1e-5], t_eval=np.linspace(0, 10, 100))
+
+        problem = pybop.FittingProblem(
+            model, parameters, dataset=dataset, signal=signal
+        )
+
+        # Test problem evaluate
+        problem_output = problem.evaluate([2e-5, 2e-5])
+
+        assert problem._model._built_model is not None
+        with pytest.raises(AssertionError):
+            np.testing.assert_allclose(
+                out["Terminal voltage [V]"].data,
+                problem_output,
+                atol=1e-5,
+            )