diff --git a/src/porepy/applications/boundary_conditions/model_boundary_conditions.py b/src/porepy/applications/boundary_conditions/model_boundary_conditions.py
index f631e6740e..4ab48ded1b 100644
--- a/src/porepy/applications/boundary_conditions/model_boundary_conditions.py
+++ b/src/porepy/applications/boundary_conditions/model_boundary_conditions.py
@@ -9,7 +9,7 @@
 import porepy as pp
 
 
-class BoundaryConditionsMassDirWestEast(pp.BoundaryConditionMixin):
+class BoundaryConditionsMassDirWestEast(pp.PorePyModel):
     """Boundary conditions for the flow problem.
 
     Dirichlet boundary conditions are defined on the west and east boundaries. Some
@@ -73,7 +73,7 @@ def bc_type_fluid_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         return pp.BoundaryCondition(sd, domain_sides.west + domain_sides.east, "dir")
 
 
-class BoundaryConditionsMassDirNorthSouth(pp.BoundaryConditionMixin):
+class BoundaryConditionsMassDirNorthSouth(pp.PorePyModel):
     """Boundary conditions for the flow problem.
 
     Dirichlet boundary conditions are defined on the north and south boundaries. Some
@@ -137,7 +137,7 @@ def bc_type_fluid_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         return pp.BoundaryCondition(sd, domain_sides.north + domain_sides.south, "dir")
 
 
-class BoundaryConditionsEnergyDirNorthSouth(pp.BoundaryConditionMixin):
+class BoundaryConditionsEnergyDirNorthSouth(pp.PorePyModel):
     """Boundary conditions for the thermal problem.
 
     Dirichlet boundary conditions are defined on the north and south boundaries. Some
@@ -183,7 +183,7 @@ def bc_type_enthalpy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         return pp.BoundaryCondition(sd, domain_sides.north + domain_sides.south, "dir")
 
 
-class BoundaryConditionsMechanicsDirNorthSouth(pp.BoundaryConditionMixin):
+class BoundaryConditionsMechanicsDirNorthSouth(pp.PorePyModel):
     """Boundary conditions for the mechanics with Dirichlet conditions on north and
     south boundaries.
 
diff --git a/src/porepy/applications/md_grids/model_geometries.py b/src/porepy/applications/md_grids/model_geometries.py
index 47b36eb5e3..9fb79f2dd1 100644
--- a/src/porepy/applications/md_grids/model_geometries.py
+++ b/src/porepy/applications/md_grids/model_geometries.py
@@ -7,7 +7,7 @@
 from . import domains, fracture_sets
 
 
-class SquareDomainOrthogonalFractures(pp.ModelGeometry):
+class SquareDomainOrthogonalFractures(pp.PorePyModel):
     """Create a mixed-dimensional grid for a square domain with up to two
     orthogonal fractures.
 
@@ -52,7 +52,7 @@ def set_domain(self) -> None:
         self._domain = domains.nd_cube_domain(2, self.domain_size)
 
 
-class CubeDomainOrthogonalFractures(pp.ModelGeometry):
+class CubeDomainOrthogonalFractures(pp.PorePyModel):
     """Create a mixed-dimensional grid for a cube domain with up to three
     orthogonal fractures.
 
@@ -78,7 +78,7 @@ def set_domain(self) -> None:
         self._domain = domains.nd_cube_domain(3, self.domain_size)
 
 
-class RectangularDomainThreeFractures(pp.ModelGeometry):
+class RectangularDomainThreeFractures(pp.PorePyModel):
     """A rectangular domain with up to three fractures.
 
     The domain is `[0, 2] x [0, 1]`.
diff --git a/src/porepy/applications/test_utils/models.py b/src/porepy/applications/test_utils/models.py
index cadecfb50e..9d8962272f 100644
--- a/src/porepy/applications/test_utils/models.py
+++ b/src/porepy/applications/test_utils/models.py
@@ -3,7 +3,7 @@
 from __future__ import annotations
 
 import inspect
-from typing import Any, Callable
+from typing import Any, Callable, cast
 
 import numpy as np
 
@@ -72,9 +72,7 @@ class Thermoporomechanics(  # type: ignore[misc]
     """Combine components needed for poromechanics simulation."""
 
 
-def model(
-    model_type: str, dim: int, num_fracs: int = 1
-) -> MassBalance | MomentumBalance | MassAndEnergyBalance | Poromechanics:
+def model(model_type: str, dim: int, num_fracs: int = 1) -> pp.PorePyModel:
     """Setup for tests."""
     # Suppress output for tests
     fracture_indices = [i for i in range(num_fracs)]
@@ -114,7 +112,7 @@ class Model(geometry, model_class):
         pass
 
     # Create an instance of the combined class
-    model = Model(params)
+    model = cast(pp.PorePyModel, Model(params))
 
     # Prepare the simulation
     # (create grids, variables, equations, discretize, etc.)
@@ -122,16 +120,10 @@ class Model(geometry, model_class):
     return model
 
 
-class RobinDirichletNeumannConditions:
+class RobinDirichletNeumannConditions(pp.PorePyModel):
     """Mixin for applying Neumann, Dirichlet and Robin conditions for a
     thermoporomechanics model."""
 
-    params: dict
-
-    domain_boundary_sides: Callable[[pp.GridLike], pp.domain.DomainSides]
-
-    nd: int
-
     def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         """Assigns pressure values on the north and south boundary."""
         p_north = self.params.get("pressure_north", 1)
@@ -259,7 +251,7 @@ def subdomains_or_interfaces_from_method_name(
     return domains
 
 
-def _add_mixin(mixin, parent):
+def _add_mixin(mixin: type, parent: type) -> type:
     """Helper method to dynamically construct a class by adding a mixin.
 
     Multiple mixins can be added by nested calls to this method.
@@ -378,7 +370,7 @@ def compare_values(
         assert np.isclose(np.sum(values_0 - values_1), 0, atol=1e-10 + rtol)
 
 
-def get_model_methods_returning_ad_operator(model_setup) -> list[str]:
+def get_model_methods_returning_ad_operator(model_setup: pp.PorePyModel) -> list[str]:
     """Get all possible testable methods to be used in test_ad_operator_methods_xx.
 
     A testable method is one that:
diff --git a/src/porepy/applications/test_utils/well_models.py b/src/porepy/applications/test_utils/well_models.py
index 6332e6a827..cc9884dbd3 100644
--- a/src/porepy/applications/test_utils/well_models.py
+++ b/src/porepy/applications/test_utils/well_models.py
@@ -1,18 +1,14 @@
 """Contains code for setting up a simple but non-trivial model with a well.
 """
 
+from typing import Literal
+
 import numpy as np
 
 import porepy as pp
 
 
-class OneVerticalWell:
-
-    domain: pp.Domain
-    """Domain for the model."""
-
-    units: pp.Units
-    """Simulation units provided by the solution strategy mixin."""
+class OneVerticalWell(pp.PorePyModel):
 
     def set_well_network(self) -> None:
         """Assign well network class."""
@@ -36,11 +32,11 @@ def meshing_arguments(self) -> dict:
 
         return mesh_sizes
 
-    def grid_type(self) -> str:
+    def grid_type(self) -> Literal["simplex"]:
         return "simplex"
 
 
-class BoundaryConditionsWellSetup(pp.BoundaryConditionMixin):
+class BoundaryConditionsWellSetup(pp.PorePyModel):
     """Boundary conditions for the well setup."""
 
     def _bc_type(self, sd: pp.Grid, well_cond: str) -> pp.BoundaryCondition:
@@ -158,6 +154,7 @@ def bc_type_fourier_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
 
 
 class WellPermeability(pp.constitutive_laws.CubicLawPermeability):
+
     def permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """Permeability [m^2].
 
diff --git a/src/porepy/examples/flow_benchmark_2d_case_4.py b/src/porepy/examples/flow_benchmark_2d_case_4.py
index 5799841b87..4b5def2fbf 100644
--- a/src/porepy/examples/flow_benchmark_2d_case_4.py
+++ b/src/porepy/examples/flow_benchmark_2d_case_4.py
@@ -20,7 +20,6 @@
 import porepy as pp
 from porepy.examples.flow_benchmark_2d_case_1 import FractureSolidConstants
 from porepy.models.constitutive_laws import DimensionDependentPermeability
-from porepy.models.protocol import PorePyModel
 
 solid_constants = FractureSolidConstants(
     residual_aperture=1e-2,  # m
@@ -30,7 +29,7 @@
 )
 
 
-class Geometry(PorePyModel):
+class Geometry(pp.PorePyModel):
     """Geometry specification."""
 
     def set_fractures(self) -> None:
@@ -43,7 +42,7 @@ def domain(self) -> pp.Domain:
         return pp.Domain({"xmax": 700, "ymax": 600})
 
 
-class BoundaryConditions(PorePyModel):
+class BoundaryConditions(pp.PorePyModel):
     """Boundary conditions for Case 4 of the 2D flow benchmark.
 
     Inflow on west (left) and prescribed pressure on east (right).
@@ -114,6 +113,6 @@ class FlowBenchmark2dCase4Model(  # type: ignore[misc]
     Geometry,
     BoundaryConditions,
     Permeability,
-    pp.fluid_mass_balance.SinglePhaseFlow,
+    pp.SinglePhaseFlow,
 ):
     """Mixer class for case 4 from the 2d flow benchmark."""
diff --git a/src/porepy/examples/flow_benchmark_3d_case_3.py b/src/porepy/examples/flow_benchmark_3d_case_3.py
index 07a8fc84a0..88de040a56 100644
--- a/src/porepy/examples/flow_benchmark_3d_case_3.py
+++ b/src/porepy/examples/flow_benchmark_3d_case_3.py
@@ -36,12 +36,9 @@
 )
 
 
-class Geometry(pp.ModelGeometry):
+class Geometry(pp.PorePyModel):
     """Define Geometry as specified in Section 5.3 of the benchmark study [1]."""
 
-    params: dict
-    """User-defined model parameters."""
-
     def set_geometry(self) -> None:
         """Create mixed-dimensional grid and fracture network."""
 
@@ -71,6 +68,7 @@ def set_geometry(self) -> None:
 
 
 class IntersectionPermeability(Permeability):
+
     def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """Constant intersection permeability.
 
@@ -90,7 +88,7 @@ def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator
         return self.isotropic_second_order_tensor(subdomains, permeability)
 
 
-class BoundaryConditions(pp.BoundaryConditionMixin):
+class BoundaryConditions(pp.PorePyModel):
     """Define inlet and outlet boundary conditions as specified by the benchmark."""
 
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
diff --git a/src/porepy/examples/mandel_biot.py b/src/porepy/examples/mandel_biot.py
index ccc3d0deaf..e390e54e38 100644
--- a/src/porepy/examples/mandel_biot.py
+++ b/src/porepy/examples/mandel_biot.py
@@ -33,13 +33,11 @@
 import scipy.optimize as opt
 
 import porepy as pp
-import porepy.models.fluid_mass_balance as mass
 import porepy.models.poromechanics as poromechanics
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.derived_models.biot import BiotPoromechanics
 from porepy.numerics.linalg.matrix_operations import sparse_array_to_row_col_data
 from porepy.utils.examples_utils import VerificationUtils
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -117,7 +115,7 @@ class MandelSaveData:
     """Current simulation time."""
 
 
-class MandelDataSaving(VerificationDataSaving):
+class MandelDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
@@ -1230,12 +1228,9 @@ def _plot_consolidation_degree(self):
 
 
 # -----> Geometry
-class MandelGeometry(pp.ModelGeometry):
+class MandelGeometry(pp.PorePyModel):
     """Class for setting up the rectangular geometry."""
 
-    params: dict
-    """Simulation model parameters."""
-
     def set_domain(self) -> None:
         """Set the domain."""
         ls = self.units.convert_units(1, "m")  # length scaling
@@ -1255,15 +1250,11 @@ def grid_type(self) -> Literal["simplex", "cartesian", "tensor_grid"]:
 
 
 # -----> Boundary conditions
-class MandelBoundaryConditionsMechanicsTimeDependent(
-    pp.momentum_balance.BoundaryConditionsMomentumBalance,
-):
+class MandelBoundaryConditionsMechanicsTimeDependent(pp.PorePyModel):
+
     exact_sol: MandelExactSolution
     """Exact solution object."""
 
-    params: dict
-    """Parameter dictionary of the verification setup."""
-
     def vertical_load(self):
         """Retrieve and scale applied force.
 
@@ -1288,8 +1279,11 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
             Vectorial boundary condition representation.
 
         """
-        # Inherit bc from parent class. This sets all bc faces as Dirichlet.
-        bc = super().bc_type_mechanics(sd=sd)
+
+        # NOTE see BC for momentum balance
+        boundary_faces = self.domain_boundary_sides(sd).all_bf
+        bc = pp.BoundaryConditionVectorial(sd, boundary_faces, "dir")
+        bc.internal_to_dirichlet(sd)
 
         # Get boundary sides, retrieve data dict, and bc object
         sides = self.domain_boundary_sides(sd)
@@ -1323,7 +1317,8 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
             the North side of the domain.
 
         """
-        bc_vals = super().bc_values_displacement(boundary_grid)
+
+        bc_vals = np.zeros((self.nd, boundary_grid.num_cells)).ravel("F")
 
         sides = self.domain_boundary_sides(boundary_grid)
         # Cells of the boundary grid are faces of the parent subdomain.
@@ -1337,7 +1332,8 @@ def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         return bc_vals
 
 
-class MandelBoundaryConditionsSinglePhaseFlow(mass.BoundaryConditionsSinglePhaseFlow):
+class MandelBoundaryConditionsSinglePhaseFlow(pp.PorePyModel):
+
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Define boundary condition types for the Darcy flux.
 
@@ -1379,21 +1375,6 @@ class MandelSolutionStrategy(poromechanics.SolutionStrategyPoromechanics):
 
     """
 
-    def __init__(self, params: dict) -> None:
-        """Constructor of the class.
-
-        Parameters:
-            params: Parameters of the verification setup.
-
-        """
-        super().__init__(params)
-
-        self.exact_sol: MandelExactSolution
-        """Exact solution object."""
-
-        self.results: list[MandelSaveData] = []
-        """List of stored results from the verification."""
-
     def set_materials(self):
         """Set material parameters.
 
diff --git a/src/porepy/examples/terzaghi_biot.py b/src/porepy/examples/terzaghi_biot.py
index 82829a3317..7313d0057e 100644
--- a/src/porepy/examples/terzaghi_biot.py
+++ b/src/porepy/examples/terzaghi_biot.py
@@ -44,13 +44,10 @@
 import numpy as np
 
 import porepy as pp
-import porepy.models.fluid_mass_balance as mass
-import porepy.models.momentum_balance as mechanics
 import porepy.models.poromechanics as poromechanics
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.derived_models.biot import BiotPoromechanics
 from porepy.utils.examples_utils import VerificationUtils
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -100,7 +97,7 @@ class TerzaghiSaveData:
     """Current simulation time."""
 
 
-class TerzaghiDataSaving(VerificationDataSaving):
+class TerzaghiDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     exact_sol: TerzaghiExactSolution
@@ -499,12 +496,9 @@ def _consolidation_degree_plot(self, color_map: mcolors.ListedColormap) -> None:
 
 
 # -----> Geometry
-class PseudoOneDimensionalColumn(pp.ModelGeometry):
+class PseudoOneDimensionalColumn(pp.PorePyModel):
     """Define geometry of the verification setup."""
 
-    params: dict
-    """Simulation model parameters."""
-
     def height(self) -> pp.number:
         """Retrieve height of the domain, in scaled [m]."""
         ls = self.units.convert_units(1, "m")  # length scaling
@@ -527,7 +521,7 @@ def meshing_arguments(self) -> dict:
 
 
 # -----> Boundary conditions
-class TerzaghiBoundaryConditionsMechanics(mechanics.BoundaryConditionsMomentumBalance):
+class TerzaghiBoundaryConditionsMechanics(pp.PorePyModel):
 
     def applied_load(self) -> pp.number:
         """Obtain vertical load in scaled [Pa]."""
@@ -545,8 +539,10 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
             the South, and rollers on the sides.
 
         """
-        # Inherit bc from parent class. This sets all bc faces as Dirichlet.
-        bc = super().bc_type_mechanics(sd=sd)
+        # Start with all faces as Dirichlet faces, analogous to base mechanics set-up.
+        boundary_faces = self.domain_boundary_sides(sd).all_bf
+        bc = pp.BoundaryConditionVectorial(sd, boundary_faces, "dir")
+        bc.internal_to_dirichlet(sd)
 
         # Get boundary sides, retrieve data dict, and bc object
         _, east, west, north, *_ = self.domain_boundary_sides(sd)
@@ -585,9 +581,8 @@ def bc_values_stress(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
         return bc_values.ravel("F")
 
 
-class TerzaghiBoundaryConditionsFlow(
-    mass.BoundaryConditionsSinglePhaseFlow,
-):
+class TerzaghiBoundaryConditionsFlow(pp.PorePyModel):
+
     def bc_type_darcy_flux(self, sd: pp.Grid) -> pp.BoundaryCondition:
         """Define boundary condition types for the Darcy flux.
 
@@ -641,21 +636,6 @@ class TerzaghiSolutionStrategy(poromechanics.SolutionStrategyPoromechanics):
 
     """
 
-    def __init__(self, params: dict) -> None:
-        """Constructor of the class.
-
-        Parameters:
-            params: Parameters of the verification setup.
-
-        """
-        super().__init__(params)
-
-        self.exact_sol: TerzaghiExactSolution
-        """Exact solution object."""
-
-        self.results: list[TerzaghiSaveData] = []
-        """List of stored results from the verification."""
-
     def set_materials(self):
         """Set material parameters.
 
diff --git a/src/porepy/models/constitutive_laws.py b/src/porepy/models/constitutive_laws.py
index 03f798c74d..102b0bbdea 100644
--- a/src/porepy/models/constitutive_laws.py
+++ b/src/porepy/models/constitutive_laws.py
@@ -3943,7 +3943,7 @@ def _mpsa_consistency(
         return consistency
 
 
-class BiotPoroMechanicsPorosity(PoroMechanicsPorosity):
+class BiotPoroMechanicsPorosity(pp.PorePyModel):
     """Porosity for poromechanical models following classical Biot's theory.
 
     The porosity is defined such that, after the chain rule is applied to the
diff --git a/src/porepy/models/derived_models/biot.py b/src/porepy/models/derived_models/biot.py
index 305e0afd47..5083d1ca99 100644
--- a/src/porepy/models/derived_models/biot.py
+++ b/src/porepy/models/derived_models/biot.py
@@ -77,11 +77,7 @@
 """
 
 import porepy as pp
-from porepy.models.poromechanics import (
-    ConstitutiveLawsPoromechanics,
-    Poromechanics,
-    SolutionStrategyPoromechanics,
-)
+from porepy.models.poromechanics import Poromechanics, SolutionStrategyPoromechanics
 
 
 class SolutionStrategyBiot(SolutionStrategyPoromechanics):
@@ -97,15 +93,9 @@ def set_materials(self):
         assert self.fluid.reference_component.compressibility == 0
 
 
-class ConstitutiveLawsBiot(
+class BiotPoromechanics(  # type: ignore[misc]
     pp.constitutive_laws.SpecificStorage,
     pp.constitutive_laws.BiotPoroMechanicsPorosity,
-    ConstitutiveLawsPoromechanics,
-): ...
-
-
-class BiotPoromechanics(  # type: ignore[misc]
-    ConstitutiveLawsBiot,
     SolutionStrategyBiot,
     Poromechanics,
 ): ...
diff --git a/src/porepy/models/protocol.py b/src/porepy/models/protocol.py
index 80fdec80eb..e6cdec136f 100644
--- a/src/porepy/models/protocol.py
+++ b/src/porepy/models/protocol.py
@@ -16,7 +16,7 @@
 """
 
 from pathlib import Path
-from typing import TYPE_CHECKING, Callable, Literal, Optional, Protocol, Sequence
+from typing import TYPE_CHECKING, Any, Callable, Literal, Optional, Protocol, Sequence
 
 import numpy as np
 import scipy.sparse as sps
@@ -76,6 +76,9 @@ def set_geometry(self) -> None:
 
             """
 
+        def set_well_network(self) -> None:
+            """Assign well network class."""
+
         def is_well(self, grid: pp.Grid | pp.MortarGrid) -> bool:
             """Check if a subdomain is a well.
 
@@ -534,6 +537,9 @@ class SolutionStrategyProtocol(Protocol):
         """Time step as an automatic differentiation scalar."""
         nonlinear_solver_statistics: pp.SolverStatistics
         """Solver statistics for the nonlinear solver."""
+        results: list[Any]
+        """A list of results collected by the data saving mixin in
+        :meth:`~porepy.viz.data_saving_model_mixin.DataSavingMixin.collect_data`."""
 
         @property
         def time_step_indices(self) -> np.ndarray:
@@ -558,6 +564,12 @@ def iterate_indices(self) -> np.ndarray:
 
             """
 
+        def prepare_simulation(self) -> None:
+            """Run at the start of simulation. Used for initialization etc."""
+
+        def after_simulation(self) -> None:
+            """Run at the end of simulation. Can be used for cleanup etc."""
+
         def _is_time_dependent(self) -> bool:
             """Specifies whether the Model problem is time-dependent.
 
@@ -699,6 +711,25 @@ def update_all_boundary_conditions(self) -> None:
 
             """
 
+        def update_boundary_condition(
+            self,
+            name: str,
+            function: Callable[[pp.BoundaryGrid], np.ndarray],
+        ) -> None:
+            """This method is the unified procedure of updating a boundary condition.
+
+            It shifts the boundary condition values in time and stores the current
+            iterate data (current time step) as the most recent previous time step data.
+            Next, it evaluates the boundary condition values for the new time step and
+            stores them in the iterate data.
+
+            Parameters:
+                name: Name of the operator defined on the boundary.
+                function: A callable that provides the boundary condition values on a
+                    given boundary grid.
+
+            """
+
     class EquationProtocol(Protocol):
         """This protocol provides declarations of methods and properties related to
         equations.
diff --git a/src/porepy/models/solution_strategy.py b/src/porepy/models/solution_strategy.py
index 41b1e93548..0db8ca44bc 100644
--- a/src/porepy/models/solution_strategy.py
+++ b/src/porepy/models/solution_strategy.py
@@ -7,7 +7,6 @@
 
 from __future__ import annotations
 
-import abc
 import logging
 import time
 import warnings
@@ -22,7 +21,7 @@
 logger = logging.getLogger(__name__)
 
 
-class SolutionStrategy(abc.ABC, pp.PorePyModel):
+class SolutionStrategy(pp.PorePyModel):
     """This is a class that specifies methods that a model must implement to
     be compatible with the linearization and time stepping methods.
 
@@ -114,6 +113,9 @@ def __init__(self, params: Optional[dict] = None):
         """Restart options. The template is provided in `SolutionStrategy.__init__`."""
         self.ad_time_step = pp.ad.Scalar(self.time_manager.dt)
         """Time step as an automatic differentiation scalar."""
+        self.results: list[Any] = []
+        """A list of results collected by the data saving mixin in
+        :meth:`~porepy.viz.data_saving_model_mixin.DataSavingMixin.collect_data`."""
 
         self.set_solver_statistics()
 
diff --git a/src/porepy/utils/porepy_types.py b/src/porepy/utils/porepy_types.py
index ae30d9ea22..1adc093417 100644
--- a/src/porepy/utils/porepy_types.py
+++ b/src/porepy/utils/porepy_types.py
@@ -5,6 +5,8 @@
 from typing import Callable, Sequence, Union
 
 import porepy as pp
+from porepy.fracs.fracture_network_2d import FractureNetwork2d
+from porepy.fracs.fracture_network_3d import FractureNetwork3d
 
 __all__ = [
     "number",
@@ -34,8 +36,8 @@
 ]
 
 fracture_network = Union[
-    "pp.fracs.fracture_network_2d.FractureNetwork2d",
-    "pp.fracs.fracture_network_3d.FractureNetwork3d",
+    "FractureNetwork2d",
+    "FractureNetwork3d",
 ]
 
 DomainFunctionType = Callable[[SubdomainsOrBoundaries], "pp.ad.Operator"]
diff --git a/src/porepy/viz/data_saving_model_mixin.py b/src/porepy/viz/data_saving_model_mixin.py
index b04a067ebd..afbcd90692 100644
--- a/src/porepy/viz/data_saving_model_mixin.py
+++ b/src/porepy/viz/data_saving_model_mixin.py
@@ -10,7 +10,7 @@
 from __future__ import annotations
 
 from pathlib import Path
-from typing import Optional, Union
+from typing import Any, Optional, Union
 
 import numpy as np
 
@@ -29,14 +29,19 @@ class DataSavingMixin(pp.PorePyModel):
     def save_data_time_step(self) -> None:
         """Export the model state at a given time step and log time.
 
-        The options for exporting times are:
-            * `None`: All time steps are exported
-            * `list`: Export if time is in the list. If the list is empty, then no
-            times are exported.
+        The options for exporting times can be given as ``params['times_to_export']``:
+
+        - ``None``: All time steps are exported.
+        - ``list``: Export if time is in the list. If the list is empty, then no
+          times are exported.
 
         In addition, save the solver statistics to file if the option is set.
 
+        Finally, :meth:`collect_data` is called and stored in :attr:`results` for
+        data collection and verification in runtime.
+
         """
+
         # Fetching the desired times to export.
         times_to_export = self.params.get("times_to_export", None)
         if times_to_export is None:
@@ -55,6 +60,40 @@ def save_data_time_step(self) -> None:
         # Save solver statistics to file.
         self.nonlinear_solver_statistics.save()
 
+        # Collecting and storing data in runtime for analysis. If default value of None
+        # is returned, nothing is stored to not burden memory.
+        if not self._is_time_dependent():  # stationary problem
+            if (
+                self.nonlinear_solver_statistics.num_iteration > 0
+            ):  # avoid saving initial condition
+                collected_data = self.collect_data()
+                if collected_data is not None:
+                    self.results.append(collected_data)
+        else:  # time-dependent problem
+            t = self.time_manager.time  # current time
+            scheduled = self.time_manager.schedule[1:]  # scheduled times except t_init
+            if any(np.isclose(t, scheduled)):
+                collected_data = self.collect_data()
+                if collected_data is not None:
+                    self.results.append(collected_data)
+
+    def collect_data(self) -> Any:
+        """Collect relevant simulation data to be stored in attr:`results`.
+
+        Override to collect data respectively. By default, this method returns None and
+        nothing is stored.
+
+        For stationary problems, this method is called in every iteration. For time
+        dependent problems, it is called after convergence of a time step which is
+        scheduled by the time manager.
+
+        Returns:
+            Any data structure relevant for future verification. By default, None.
+            If it is not None, it is stored in :attr:`results`.
+
+        """
+        return None
+
     def write_pvd_and_vtu(self) -> None:
         """Helper function for writing the .vtu and .pvd files and time information."""
         self.exporter.write_vtu(self.data_to_export(), time_dependent=True)
@@ -237,29 +276,3 @@ def load_data_from_pvd(
         self.time_manager.load_time_information(times_file)
         self.time_manager.set_time_and_dt_from_exported_steps(time_index)
         self.exporter._time_step_counter = time_index
-
-
-class VerificationDataSaving(DataSavingMixin):
-    """Class to store relevant data for a generic verification setup."""
-
-    results: list
-    """List of objects containing the results of the verification."""
-
-    def save_data_time_step(self) -> None:
-        """Save data to the `results` list."""
-        if not self._is_time_dependent():  # stationary problem
-            if (
-                self.nonlinear_solver_statistics.num_iteration > 0
-            ):  # avoid saving initial condition
-                collected_data = self.collect_data()
-                self.results.append(collected_data)
-        else:  # time-dependent problem
-            t = self.time_manager.time  # current time
-            scheduled = self.time_manager.schedule[1:]  # scheduled times except t_init
-            if any(np.isclose(t, scheduled)):
-                collected_data = self.collect_data()
-                self.results.append(collected_data)
-
-    def collect_data(self):
-        """Collect relevant data for the verification setup."""
-        raise NotImplementedError()
diff --git a/tests/applications/test_convergence_analysis.py b/tests/applications/test_convergence_analysis.py
index c511532be0..7d30e4c4ef 100644
--- a/tests/applications/test_convergence_analysis.py
+++ b/tests/applications/test_convergence_analysis.py
@@ -26,7 +26,6 @@
 )
 from porepy.models.fluid_mass_balance import SinglePhaseFlow
 from porepy.utils.txt_io import read_data_from_txt
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 
 # -----> Fixtures that are required on a module level.
@@ -461,7 +460,7 @@ class StationaryModelSaveData:
         error_var_0: float  # error associated with variable 0
         error_var_1: float  # error associated with variable 1
 
-    class StationaryModelDataSaving(VerificationDataSaving):
+    class StationaryModelDataSaving(pp.PorePyModel):
         """Class that collects and store data."""
 
         def collect_data(self) -> StationaryModelSaveData:
@@ -485,9 +484,7 @@ def collect_data(self) -> StationaryModelSaveData:
     class StationaryModelSolutionStrategy(pp.SolutionStrategy):
         """Solution strategy for the stationary flow model."""
 
-        def __init__(self, params: dict):
-            super().__init__(params)
-            self.results: list[StationaryModelSaveData] = []
+        results: list[StationaryModelSaveData]
 
         def _is_nonlinear_problem(self) -> bool:
             """Whether the model is non-linear."""
@@ -528,7 +525,7 @@ class TimeDependentModelSaveData:
         error_var_0: float  # error associated with variable 0
         error_var_1: float  # error associated with variable 1
 
-    class TimeDependentModelDataSaving(VerificationDataSaving):
+    class TimeDependentModelDataSaving(pp.PorePyModel):
         """Class that collects and store data."""
 
         def collect_data(self) -> TimeDependentModelSaveData:
@@ -552,9 +549,7 @@ def collect_data(self) -> TimeDependentModelSaveData:
     class TimeDependentModelSolutionStrategy(pp.SolutionStrategy):
         """Solution strategy for the time-dependent flow model."""
 
-        def __init__(self, params: dict):
-            super().__init__(params)
-            self.results: list[TimeDependentModelSaveData] = []
+        results: list[TimeDependentModelSaveData]
 
         def _is_nonlinear_problem(self) -> bool:
             """Whether the problem is non-linear."""
diff --git a/tests/functional/setups/manu_flow_comp_2d_frac.py b/tests/functional/setups/manu_flow_comp_2d_frac.py
index ae0e14ecd8..a6511722c0 100644
--- a/tests/functional/setups/manu_flow_comp_2d_frac.py
+++ b/tests/functional/setups/manu_flow_comp_2d_frac.py
@@ -34,7 +34,6 @@
 
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 from tests.functional.setups.manu_flow_incomp_frac_2d import (
     ManuIncompSaveData,
     ManuIncompUtils,
@@ -75,7 +74,7 @@ class ManuCompSaveData(ManuIncompSaveData):
     """Current simulation time."""
 
 
-class ManuCompDataSaving(VerificationDataSaving):
+class ManuCompDataSaving(pp.PorePyModel):
     """Mixin class to store relevant data."""
 
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
@@ -703,12 +702,6 @@ def __init__(self, params: dict):
         """Constructor of the class."""
         super().__init__(params)
 
-        self.exact_sol: ManuCompExactSolution2d
-        """Exact solution object."""
-
-        self.results: list[ManuCompSaveData] = []
-        """Object that stores exact and approximated solutions and L2 errors."""
-
         self.subdomain_darcy_flux_variable: str = "darcy_flux"
         """Keyword to access the subdomain Darcy fluxes."""
 
diff --git a/tests/functional/setups/manu_flow_incomp_frac_2d.py b/tests/functional/setups/manu_flow_incomp_frac_2d.py
index e2edb1c694..9ae5d92de8 100644
--- a/tests/functional/setups/manu_flow_incomp_frac_2d.py
+++ b/tests/functional/setups/manu_flow_incomp_frac_2d.py
@@ -25,7 +25,6 @@
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.applications.md_grids.domains import nd_cube_domain
 from porepy.utils.examples_utils import VerificationUtils
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -97,7 +96,7 @@ class ManuIncompSaveData:
     """Exact pressure in the matrix."""
 
 
-class ManuIncompDataSaving(VerificationDataSaving):
+class ManuIncompDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
@@ -755,17 +754,6 @@ class ManuIncompSolutionStrategy2d(
     results: list[ManuIncompSaveData]
     """List of SaveData objects."""
 
-    def __init__(self, params: dict):
-        """Constructor for the class."""
-
-        super().__init__(params)
-
-        self.exact_sol: ManuIncompExactSolution2d
-        """Exact solution object."""
-
-        self.results: list[ManuIncompSaveData] = []
-        """Results object that stores exact and approximated solutions and errors."""
-
     def set_materials(self):
         """Set material constants for the verification setup."""
         super().set_materials()
diff --git a/tests/functional/setups/manu_poromech_nofrac_2d.py b/tests/functional/setups/manu_poromech_nofrac_2d.py
index 193727b931..0cebdf2daa 100644
--- a/tests/functional/setups/manu_poromech_nofrac_2d.py
+++ b/tests/functional/setups/manu_poromech_nofrac_2d.py
@@ -60,7 +60,6 @@
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.applications.md_grids.domains import nd_cube_domain
 from porepy.utils.examples_utils import VerificationUtils
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -112,7 +111,7 @@ class ManuPoroMechSaveData:
     """Current simulation time."""
 
 
-class ManuPoroMechDataSaving(VerificationDataSaving):
+class ManuPoroMechDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
@@ -717,12 +716,6 @@ def __init__(self, params: dict):
         """Constructor for the class."""
         super().__init__(params)
 
-        self.exact_sol: ManuPoroMechExactSolution2d
-        """Exact solution object."""
-
-        self.results: list[ManuPoroMechSaveData] = []
-        """Results object that stores exact and approximated solutions and errors."""
-
         self.flux_variable: str = "darcy_flux"
         """Keyword to access the Darcy fluxes."""
 
diff --git a/tests/functional/setups/manu_thermoporomech_nofrac_2d.py b/tests/functional/setups/manu_thermoporomech_nofrac_2d.py
index 892839a003..54c1293700 100644
--- a/tests/functional/setups/manu_thermoporomech_nofrac_2d.py
+++ b/tests/functional/setups/manu_thermoporomech_nofrac_2d.py
@@ -47,7 +47,6 @@
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.applications.md_grids.domains import nd_cube_domain
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 
 # PorePy typings
 number = pp.number
@@ -116,7 +115,7 @@ class ManuThermoPoroMechSaveData:
     time: number
 
 
-class ManuThermoPoroMechDataSaving(VerificationDataSaving):
+class ManuThermoPoroMechDataSaving(pp.PorePyModel):
     """Mixin class to save relevant data."""
 
     exact_sol: ManuThermoPoroMechExactSolution2d
@@ -138,6 +137,12 @@ class ManuThermoPoroMechDataSaving(VerificationDataSaving):
 
     """
 
+    energy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Method that returns the energy fluxes in the form of an Ad operator. Usually
+    provided by the mixin class
+    :class:`porepy.models.energy_balance.EnergyBalanceEquations`.
+
+    """
     darcy_flux: Callable[[list[pp.Grid]], pp.ad.Operator]
     """Method that returns the Darcy fluxes in the form of an Ad operator. Usually
     provided by the mixin class :class:`porepy.models.constitutive_laws.DarcysLaw`.
@@ -895,7 +900,7 @@ def energy_source(self, sd: pp.Grid, time: float) -> np.ndarray:
 
 
 # -----> Geometry
-class UnitSquareGrid(pp.ModelGeometry):
+class UnitSquareGrid(pp.PorePyModel):
     """Class for setting up the geometry of the unit square domain.
 
     The domain may be assigned different material parameters in the region x > 0.5 and y
@@ -917,9 +922,6 @@ class UnitSquareGrid(pp.ModelGeometry):
 
     """
 
-    params: dict
-    """Simulation model parameters."""
-
     def set_geometry(self) -> None:
         super().set_geometry()
 
@@ -1034,16 +1036,13 @@ class ManuThermoPoroMechSolutionStrategy2d(
 ):
     """Solution strategy for the verification setup."""
 
+    exact_sol: ManuThermoPoroMechExactSolution2d
+    """Exact solution object."""
+
     def __init__(self, params: dict):
         """Constructor for the class."""
         super().__init__(params)
 
-        self.exact_sol: ManuThermoPoroMechExactSolution2d
-        """Exact solution object."""
-
-        self.results: list[ManuThermoPoroMechSaveData] = []
-        """Results object that stores exact and approximated solutions and errors."""
-
         self.flux_variable: str = "darcy_flux"
         """Keyword to access the Darcy fluxes."""
 
diff --git a/tests/functional/setups/manu_thermoporomech_nofrac_3d.py b/tests/functional/setups/manu_thermoporomech_nofrac_3d.py
index 4b988c2cf5..9a02da7a06 100644
--- a/tests/functional/setups/manu_thermoporomech_nofrac_3d.py
+++ b/tests/functional/setups/manu_thermoporomech_nofrac_3d.py
@@ -791,7 +791,7 @@ def energy_source(self, sd: pp.Grid, time: float) -> np.ndarray:
 
 
 # -----> Geometry
-class UnitCubeGrid(pp.ModelGeometry):
+class UnitCubeGrid(pp.PorePyModel):
     """Class for setting up the geometry of the unit cube domain.
 
     The domain may be assigned different material parameters in the region x > 0.5, y >
@@ -812,9 +812,6 @@ class UnitCubeGrid(pp.ModelGeometry):
 
     """
 
-    params: dict
-    """Simulation model parameters."""
-
     def set_geometry(self) -> None:
         super().set_geometry()
 
@@ -876,18 +873,15 @@ class ManuThermoPoroMechSolutionStrategy3d(
 ):
     """Solution strategy for the verification setup."""
 
+    exact_sol: ManuThermoPoroMechExactSolution3d
+    """Exact solution object."""
+
     def __init__(self, params: dict):
         """Constructor for the class."""
         super().__init__(params)
 
-        self.exact_sol: ManuThermoPoroMechExactSolution3d
-        """Exact solution object."""
-
         self.stress_variable: str = "thermoporoelastic_force"
-        """Keyword to access the thermoporoelastic force."""
-        
-        self.results: list[ManuThermoPoroMechSaveData] = []
-        """Results object that stores exact and approximated solutions and errors."""        
+        """Keyword to access the thermoporoelastic force."""      
 
     def set_materials(self):
         """Set material parameters."""
diff --git a/tests/functional/test_benchmark_2d_case_3.py b/tests/functional/test_benchmark_2d_case_3.py
index 55b1d0224c..54c625a931 100644
--- a/tests/functional/test_benchmark_2d_case_3.py
+++ b/tests/functional/test_benchmark_2d_case_3.py
@@ -71,6 +71,7 @@ def model(
         "grid_type": "simplex",
         "meshing_arguments": {"cell_size": 0.1},
         "flux_discretization": flux_discretization,
+        "times_to_export": [],  # Suppress output for tests
     }
     if case == "a":
         model = Model3aWithEffectivePermeability(model_params)
diff --git a/tests/functional/test_benchmark_3d_case_3.py b/tests/functional/test_benchmark_3d_case_3.py
index 4e7dad4586..60a95e51c2 100644
--- a/tests/functional/test_benchmark_3d_case_3.py
+++ b/tests/functional/test_benchmark_3d_case_3.py
@@ -53,6 +53,7 @@ def model(
     model_params = {
         "material_constants": {"solid": solid_constants},
         "flux_discretization": flux_discretization,
+        "times_to_export": [],  # Suppress output for tests
     }
     model = ModelWithEffectivePermeability(model_params)
     pp.run_time_dependent_model(model)
diff --git a/tests/functional/test_mandel.py b/tests/functional/test_mandel.py
index db612d1ded..d5566c93d6 100644
--- a/tests/functional/test_mandel.py
+++ b/tests/functional/test_mandel.py
@@ -41,6 +41,7 @@ def results() -> list[MandelSaveData]:
     model_params = {
         "material_constants": material_constants,
         "time_manager": time_manager,
+        "times_to_export": [],  # Suppress output for tests
     }
     setup = MandelSetup(model_params)
     pp.run_time_dependent_model(setup)
@@ -78,7 +79,9 @@ def results() -> list[MandelSaveData]:
 
 
 @pytest.mark.parametrize("time_index", [0, 1])
-def test_error_primary_and_secondary_variables(time_index: int, results):
+def test_error_primary_and_secondary_variables(
+    time_index: int, results: list[MandelSaveData]
+):
     """Checks error for pressure, displacement, flux, force, and consolidation degree.
 
     Physical parameters used in this test have been adapted from [1].
@@ -153,6 +156,7 @@ def test_scaled_vs_unscaled_systems():
     model_params_unscaled = {
         "material_constants": material_constants_unscaled,
         "time_manager": time_manager_unscaled,
+        "times_to_export": [],  # Suppress output for tests
     }
     model_unscaled = MandelSetup(params=model_params_unscaled)
     pp.run_time_dependent_model(model_unscaled)
@@ -169,6 +173,7 @@ def test_scaled_vs_unscaled_systems():
         "material_constants": material_constants_scaled,
         "time_manager": time_manager_scaled,
         "units": units,
+        "times_to_export": [],  # Suppress output for tests
     }
     scaled_model = MandelSetup(params=model_params_scaled)
     pp.run_time_dependent_model(model=scaled_model)
diff --git a/tests/functional/test_manu_flow_comp_frac.py b/tests/functional/test_manu_flow_comp_frac.py
index 8b5858f5d9..844c990c67 100644
--- a/tests/functional/test_manu_flow_comp_frac.py
+++ b/tests/functional/test_manu_flow_comp_frac.py
@@ -109,6 +109,7 @@ def actual_l2_errors(
         "reference_variable_values": reference_values,
         "meshing_arguments": {"cell_size": 0.125},
         "time_manager": pp.TimeManager([0, 0.5, 1.0], 0.5, True),
+        "times_to_export": [],  # Suppress output for tests
     }
 
     # Retrieve actual L2-relative errors.
@@ -116,7 +117,7 @@ def actual_l2_errors(
     # Loop through models, i.e., 2d and 3d.
     for model in [ManuCompFlowSetup2d, ManuCompFlowSetup3d]:
         # Make deep copy of params to avoid nasty bugs.
-        setup = model(deepcopy(model_params))
+        setup: pp.PorePyModel = model(deepcopy(model_params))
         pp.run_time_dependent_model(setup, {})
         errors_setup: list[dict[str, float]] = []
         # Loop through results, i.e., results for each scheduled time.
@@ -279,6 +280,7 @@ def actual_ooc(
                     "material_constants": material_constants,
                     "reference_variable_values": reference_values,
                     "meshing_arguments": {"cell_size": 0.125},
+                    "times_to_export": [],  # Suppress output for tests
                 }
                 # Use 4 levels of refinement for 2d and 3 levels for 3d
                 if model_idx == 0:
diff --git a/tests/functional/test_manu_flow_incomp_frac.py b/tests/functional/test_manu_flow_incomp_frac.py
index a9282a9618..bf3322900e 100644
--- a/tests/functional/test_manu_flow_incomp_frac.py
+++ b/tests/functional/test_manu_flow_incomp_frac.py
@@ -82,6 +82,7 @@ def actual_l2_errors(material_constants: dict) -> list[dict[str, float]]:
         "grid_type": "cartesian",
         "material_constants": material_constants,
         "meshing_arguments": {"cell_size": 0.125},
+        "times_to_export": [],  # Suppress output for tests
     }
 
     # Retrieve actual L2-relative errors
@@ -89,7 +90,9 @@ def actual_l2_errors(material_constants: dict) -> list[dict[str, float]]:
     # Loop through models, i.e., 2d and 3d
     for model in [ManuIncompFlowSetup2d, ManuIncompFlowSetup3d]:
         # Make deep copy of params to avoid nasty bugs.
-        setup = model(deepcopy(model_params))
+        setup: ManuIncompFlowSetup2d | ManuIncompFlowSetup3d = model(
+            deepcopy(model_params)
+        )
         pp.run_time_dependent_model(setup)
         errors.append(
             {
@@ -217,6 +220,7 @@ def actual_ooc(material_constants: dict) -> list[list[dict[str, float]]]:
                     "grid_type": grid_type,
                     "material_constants": material_constants,
                     "meshing_arguments": {"cell_size": 0.125},
+                    "times_to_export": [],  # Suppress output for tests
                 }
                 # Use 4 levels of refinement for 2d and 3 levels for 3d
                 if model_idx == 0:
diff --git a/tests/functional/test_manu_poromech_nofrac.py b/tests/functional/test_manu_poromech_nofrac.py
index a1f3dbf7b0..6db8fa0b37 100644
--- a/tests/functional/test_manu_poromech_nofrac.py
+++ b/tests/functional/test_manu_poromech_nofrac.py
@@ -96,6 +96,7 @@ def actual_l2_errors(material_constants: dict) -> list[list[dict[str, float]]]:
         "meshing_arguments": {"cell_size": 0.25},
         "manufactured_solution": "nordbotten_2016",
         "time_manager": pp.TimeManager([0, 0.5, 1.0], 0.5, True),
+        "times_to_export": [],  # Suppress output for tests
     }
 
     # Retrieve actual L2-relative errors.
@@ -103,7 +104,7 @@ def actual_l2_errors(material_constants: dict) -> list[list[dict[str, float]]]:
     # Loop through models, i.e., 2d and 3d.
     for model in [ManuPoroMechSetup2d, ManuPoroMechSetup3d]:
         # Make deep copy of params to avoid nasty bugs.
-        setup = model(deepcopy(model_params))
+        setup: pp.PorePyModel = model(deepcopy(model_params))
         pp.run_time_dependent_model(setup)
         errors_setup: list[dict[str, float]] = []
         # Loop through results, i.e., results for each scheduled time.
@@ -251,6 +252,7 @@ def actual_ooc(material_constants: dict) -> list[list[dict[str, float]]]:
                     "grid_type": grid_type,
                     "material_constants": material_constants,
                     "meshing_arguments": {"cell_size": 0.25},
+                    "times_to_export": [],  # Suppress output for tests
                 }
                 # Use 4 levels of refinement for 2d and 3 levels for 3d.
                 if model_idx == 0:
diff --git a/tests/functional/test_manu_thermoporomech_nofrac.py b/tests/functional/test_manu_thermoporomech_nofrac.py
index da4f4bc0f4..6e89f9e043 100644
--- a/tests/functional/test_manu_thermoporomech_nofrac.py
+++ b/tests/functional/test_manu_thermoporomech_nofrac.py
@@ -82,6 +82,7 @@ def actual_l2_errors(material_constants) -> list[list[dict[str, float]]]:
         "meshing_arguments": {"cell_size": 0.25},
         "time_manager": pp.TimeManager([0, 0.5, 1.0], 0.5, True),
         "heterogeneity": 10.0,
+        "times_to_export": [],  # Suppress output for tests
     }
 
     # Retrieve actual L2-relative errors.
@@ -89,7 +90,7 @@ def actual_l2_errors(material_constants) -> list[list[dict[str, float]]]:
     # Loop through models, i.e., 2d and 3d.
     for model in [ManuThermoPoroMechSetup2d, ManuThermoPoroMechSetup3d]:
         # Make deep copy of params to avoid nasty bugs.
-        setup = model(deepcopy(model_params))
+        setup: pp.PorePyModel = model(deepcopy(model_params))
         pp.run_time_dependent_model(setup, {})
         errors_setup: list[dict[str, float]] = []
         # Loop through results, i.e., results for each scheduled time.
@@ -253,6 +254,7 @@ def actual_ooc(material_constants: dict) -> list[list[dict[str, float]]]:
                     "meshing_arguments": {"cell_size": 0.25},
                     "perturbation": 0.3,
                     "heterogeneity": 10.0,
+                    "times_to_export": [],  # Suppress output for tests
                 }
                 # Use 4 levels of refinement for 2d and 3 levels for 3d.
                 if model_idx == 0:
diff --git a/tests/functional/test_terzaghi.py b/tests/functional/test_terzaghi.py
index e2b06402d4..fe909e9777 100644
--- a/tests/functional/test_terzaghi.py
+++ b/tests/functional/test_terzaghi.py
@@ -56,6 +56,7 @@ def test_biot_equal_to_incompressible_poromechanics():
             "fluid": pp.FluidComponent(**terzaghi_fluid_constants),
         },
         "num_cells": 10,
+        "times_to_export": [],  # Suppress output for tests
     }
     setup_poromech = TerzaghiSetupPoromechanics(model_params_poromech)
     pp.run_time_dependent_model(model=setup_poromech)
@@ -69,6 +70,7 @@ def test_biot_equal_to_incompressible_poromechanics():
             "fluid": pp.FluidComponent(**terzaghi_fluid_constants),
         },
         "num_cells": 10,
+        "times_to_export": [],  # Suppress output for tests
     }
     setup_biot = TerzaghiSetup(model_params_biot)
     pp.run_time_dependent_model(model=setup_biot)
@@ -114,6 +116,7 @@ def test_pressure_and_consolidation_degree_errors():
         },
         "time_manager": pp.TimeManager([0, 0.15, 0.3], 0.15, True),
         "num_cells": 10,
+        "times_to_export": [],  # Suppress output for tests
     }
     setup = TerzaghiSetup(model_params)
     pp.run_time_dependent_model(setup)
@@ -144,7 +147,10 @@ def test_scaled_vs_unscaled_systems():
         "fluid": pp.FluidComponent(**terzaghi_fluid_constants),
         "solid": pp.SolidConstants(**terzaghi_solid_constants),
     }
-    model_params_unscaled = {"material_constants": material_constants_unscaled}
+    model_params_unscaled = {
+        "material_constants": material_constants_unscaled,
+        "times_to_export": [],  # Suppress output for tests
+    }
     unscaled = TerzaghiSetup(params=model_params_unscaled)
     pp.run_time_dependent_model(model=unscaled)
 
@@ -155,7 +161,11 @@ def test_scaled_vs_unscaled_systems():
     }
     scaling = {"m": 0.001, "kg": 0.001}  # length in millimeters and mass in grams
     units = pp.Units(**scaling)
-    model_params_scaled = {"material_constants": material_constants_scaled, "units": units}
+    model_params_scaled = {
+        "material_constants": material_constants_scaled,
+        "units": units,
+        "times_to_export": [],  # Suppress output for tests
+    }
     scaled = TerzaghiSetup(params=model_params_scaled)
     pp.run_time_dependent_model(model=scaled)
 
diff --git a/tests/models/test_boundary_condition.py b/tests/models/test_boundary_condition.py
index ce11a43295..8da495fec1 100644
--- a/tests/models/test_boundary_condition.py
+++ b/tests/models/test_boundary_condition.py
@@ -15,7 +15,7 @@
 from porepy.models.momentum_balance import MomentumBalance
 
 
-class CustomBoundaryCondition(pp.BoundaryConditionMixin):
+class CustomBoundaryCondition(pp.PorePyModel):
     """We define a custom dummy boundary condition.
 
     Neumann values are explicitly set, they are time dependent.
@@ -76,7 +76,11 @@ def test_boundary_condition_mixin(t_end: int):
     3) Previous timestep values are set correctly for the time dependent Neumann.
 
     """
-    setup = MassBalance()
+    setup = MassBalance(
+        {
+            "times_to_export": [],  # Suppress output for tests
+        }
+    )
     setup.time_manager.dt = 1
     setup.time_manager.time_final = t_end
     pp.run_time_dependent_model(setup)
@@ -385,6 +389,7 @@ def run_model(balance_class, alpha):
         "times_to_export": [],
         "fracture_indices": [],
         "meshing_arguments": {"cell_size": 0.5},
+        "times_to_export": [],  # Suppress output for tests
     }
 
     params["alpha"] = alpha
diff --git a/tests/models/test_constitutive_laws.py b/tests/models/test_constitutive_laws.py
index 0f23cdd372..3573d5db42 100644
--- a/tests/models/test_constitutive_laws.py
+++ b/tests/models/test_constitutive_laws.py
@@ -348,6 +348,7 @@ def test_evaluated_values(
     params = {
         "material_constants": {"solid": solid, "fluid": fluid},
         "fracture_indices": [0, 1],
+        "times_to_export": [],  # Suppress output for tests
     }
 
     setup = model(params)
@@ -419,6 +420,7 @@ def test_perturbation_from_reference(model: type[models.MassAndEnergyBalance], q
         "reference_variable_values": pp.ReferenceVariableValues(
             pressure=1, temperature=2
         ),
+        "times_to_export": [],
     }
 
     setup = model(params)
@@ -466,7 +468,11 @@ def test_dimension_reduction_values(
     """
     # Assign non-trivial values to the parameters to avoid masking errors.
     solid = pp.SolidConstants(residual_aperture=0.02)
-    params = {"material_constants": {"solid": solid}, "num_fracs": 3}
+    params = {
+        "material_constants": {"solid": solid},
+        "num_fracs": 3,
+        "times_to_export": [],
+    }
     if geometry is models.RectangularDomainThreeFractures:
         params["fracture_indices"] = [0, 1]
 
@@ -635,7 +641,9 @@ def test_derivatives_darcy_flux_potential_trace(base_discr: str):
     """
 
     # Set up and discretize model
-    model = PoromechanicalTestDiffTpfa({"darcy_flux_discretization": base_discr})
+    model = PoromechanicalTestDiffTpfa(
+        {"darcy_flux_discretization": base_discr, "times_to_export": []}
+    )
     model.prepare_simulation()
     model.discretize()
 
diff --git a/tests/models/test_energy_balance.py b/tests/models/test_energy_balance.py
index 78fa264828..a71ccbba5e 100644
--- a/tests/models/test_energy_balance.py
+++ b/tests/models/test_energy_balance.py
@@ -231,7 +231,7 @@ def bc_values_pressure(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:
     # Non-unitary time step needed for convergence
     dt = 1e5
     model_params = {
-        "times_to_export": [],
+        "times_to_export": [],  # Suppress output for tests
         "fracture_indices": [0, 1],
         "cartesian": True,
         "material_constants": {"solid": solid, "fluid": fluid, "numerical": numerical},
@@ -308,6 +308,7 @@ def test_energy_conservation():
         "fracture_indices": [2],
         "time_manager": pp.TimeManager(schedule=[0, dt], dt_init=dt, constant_dt=True),
         "grid_type": "cartesian",
+        "times_to_export": [],  # Suppress output for tests
     }
 
     setup = MassAndEnergyWellModel(model_params)
diff --git a/tests/models/test_fluid_mass_balance.py b/tests/models/test_fluid_mass_balance.py
index 09d508d635..13ad862937 100644
--- a/tests/models/test_fluid_mass_balance.py
+++ b/tests/models/test_fluid_mass_balance.py
@@ -74,6 +74,7 @@ class Model(SquareDomainOrthogonalFractures, SinglePhaseFlow):
         "fracture_indices": [0, 1],
         "grid_type": "cartesian",
         "meshing_arguments": {"cell_size_x": 0.5, "cell_size_y": 0.5},
+        "times_to_export": [],
     }
 
     # Instantiate the model setup
@@ -550,6 +551,7 @@ def test_well_incompressible_pressure_values():
         },
         # Use only the horizontal fracture of OrthogonalFractures3d
         "fracture_indices": [2],
+        "times_to_export": [],
     }
 
     setup = WellModel(params)
@@ -670,6 +672,7 @@ def model_setup_gravity(
         "fracture_indices": [-1],  # Constant y and z coordinates in 2d and 3d, resp.
         "meshing_arguments": {"cell_size": 0.5},
         "darcy_flux_discretization": discretization_method,
+        "times_to_export": [],
     }
     params.update(model_params)
     params["material_constants"] = {
@@ -893,6 +896,7 @@ def test_no_flow_neumann(
         params = {
             "meshing_arguments": {"cell_size": 1 / 3, "cell_size_y": 1 / 2},
             "grid_type": grid_type,
+            "times_to_export": [],
         }
         self.model = model_setup_gravity(
             dimension=2,
@@ -931,6 +935,7 @@ def test_no_flow_rotate_gravity(self, discretization, num_nodes_mortar, angle):
         params = {
             "meshing_arguments": {"cell_size": 1.0, "cell_size_y": 1 / 2},
             "grid_type": "cartesian",
+            "times_to_export": [],
         }
         num_nodes_1d = 2
         self.model = model_setup_gravity(
@@ -965,6 +970,7 @@ def test_no_flow_dirichlet(
         params = {
             "meshing_arguments": {"cell_size": 1 / 3, "cell_size_y": 1 / 2},
             "grid_type": grid_type,
+            "times_to_export": [],
         }
 
         self.model = model_setup_gravity(
@@ -993,6 +999,7 @@ def test_inflow_top(
         params = {
             "meshing_arguments": {"cell_size": 1 / 2},
             "grid_type": grid_type,
+            "times_to_export": [],
         }
         a = 1e-2
         self.model = model_setup_gravity(
diff --git a/tests/models/test_geometry.py b/tests/models/test_geometry.py
index 4014793ef3..5e5015eacc 100644
--- a/tests/models/test_geometry.py
+++ b/tests/models/test_geometry.py
@@ -27,8 +27,12 @@
 
 
 # List of geometry classes to test.
+# Turn mixins of specific grids into proper model geometries.
 geometry_list: list[pp.ModelGeometry] = [
-    porepy.applications.md_grids.model_geometries.RectangularDomainThreeFractures,
+    models._add_mixin(
+        porepy.applications.md_grids.model_geometries.RectangularDomainThreeFractures,
+        pp.ModelGeometry,
+    ),
     models._add_mixin(
         porepy.applications.md_grids.model_geometries.OrthogonalFractures3d,
         pp.ModelGeometry,
diff --git a/tests/models/test_poromechanics.py b/tests/models/test_poromechanics.py
index 2689b0f262..364ab9f2d8 100644
--- a/tests/models/test_poromechanics.py
+++ b/tests/models/test_poromechanics.py
@@ -317,7 +317,7 @@ def test_poromechanics_model_no_modification():
 
     Failure of this test would signify rather fundamental problems in the model.
     """
-    mod = pp.Poromechanics({})
+    mod = pp.Poromechanics({"times_to_export": []})
     pp.run_stationary_model(mod, {})
 
 
@@ -330,6 +330,7 @@ def test_without_fracture(biot_coefficient):
         "material_constants": {"fluid": fluid, "solid": solid},
         "u_north": [0.0, 0.001],
         "cartesian": True,
+        "times_to_export": [],
     }
     m = TailoredPoromechanics(params)
     pp.run_time_dependent_model(m)
@@ -539,6 +540,7 @@ def test_poromechanics_well():
     model_params = {
         "fracture_indices": [2],
         "well_flux": -1e-2,
+        "times_to_export": [],
     }
     setup = PoromechanicsWell(model_params)
     pp.run_time_dependent_model(setup)
diff --git a/tests/models/test_solution_strategy.py b/tests/models/test_solution_strategy.py
index c16f724786..1dc89f6323 100644
--- a/tests/models/test_solution_strategy.py
+++ b/tests/models/test_solution_strategy.py
@@ -264,6 +264,7 @@ def test_targeted_rediscretization(model_class):
         "cartesian": True,
         # Make flow problem non-linear:
         "material_constants": {"fluid": pp.FluidComponent(compressibility=1.0)},
+        "times_to_export": [],
     }
     # Finalize the model class by adding the rediscretization mixin.
     rediscretization_model_class = models._add_mixin(RediscretizationTest, model_class)
diff --git a/tests/models/test_thermoporomechanics.py b/tests/models/test_thermoporomechanics.py
index 344b8a9d3e..7ef751cf17 100644
--- a/tests/models/test_thermoporomechanics.py
+++ b/tests/models/test_thermoporomechanics.py
@@ -178,7 +178,7 @@ def test_thermoporomechanics_model_no_modification():
     Failure of this test would signify rather fundamental problems in the model.
 
     """
-    mod = pp.Thermoporomechanics({})
+    mod = pp.Thermoporomechanics({"times_to_export": []})
     pp.run_stationary_model(mod, {})
 
 
@@ -304,6 +304,7 @@ def set_domain(self) -> None:
         "fourier_flux_east": -2e-2,
         "mechanical_stress_west": 3e-2,
         "mechanical_stress_east": -3e-2,
+        "times_to_export": [],
     }
 
     model = TailoredPoromechanicsRobin(model_params)
@@ -463,6 +464,7 @@ def test_thermoporomechanics_well():
     model_params = {
         "fracture_indices": [2],
         "well_flux": -1e-2,
+        "times_to_export": [],
     }
     setup = ThermoporomechanicsWell(model_params)
     pp.run_time_dependent_model(setup)
diff --git a/tests/numerics/fracture_deformation/test_fracture_propagation.py b/tests/numerics/fracture_deformation/test_fracture_propagation.py
index 976dbe2121..5356b037c8 100644
--- a/tests/numerics/fracture_deformation/test_fracture_propagation.py
+++ b/tests/numerics/fracture_deformation/test_fracture_propagation.py
@@ -166,7 +166,7 @@ def test_pick_propagation_face_conforming_propagation(case):
     data_primary[pp.TIME_STEP_SOLUTIONS] = {}
 
     # Propagation model; assign this some necessary fields.
-    model = pp.ConformingFracturePropagation({})
+    model = pp.ConformingFracturePropagation({"times_to_export": []})
     model.mechanics_parameter_key = mech_key
     model.nd = mdg.dim_max()
 
diff --git a/tests/numerics/fv/test_tpfa.py b/tests/numerics/fv/test_tpfa.py
index 4f46ee29ae..db384fcca2 100644
--- a/tests/numerics/fv/test_tpfa.py
+++ b/tests/numerics/fv/test_tpfa.py
@@ -333,6 +333,7 @@ def test_transmissibility_calculation(vector_source: bool, base_discr: str):
     model_params = {
         "darcy_flux_discretization": base_discr,
         "vector_source": vector_source_array,
+        "times_to_export": [],
     }
 
     model = UnitTestAdTpfaFlux(model_params)
@@ -671,7 +672,11 @@ def test_diff_tpfa_on_grid_with_all_dimensions(base_discr: str, grid_type: str):
 
     """
     model = DiffTpfaGridsOfAllDimensions(
-        {"darcy_flux_discretization": base_discr, "grid_type": grid_type}
+        {
+            "darcy_flux_discretization": base_discr,
+            "grid_type": grid_type,
+            "times_to_export": []
+        }
     )
     model.prepare_simulation()
 
@@ -762,6 +767,7 @@ def test_diff_tpfa_and_standard_tpfa_give_same_linear_system(base_discr: str):
     params = {
         "darcy_flux_discretization": base_discr,
         "fourier_flux_discretization": base_discr,
+        "times_to_export": [],
     }
     model_without_diff = WithoutDiffTpfa(params.copy())
     model_with_diff = WithDiffTpfa(params)
@@ -838,7 +844,9 @@ def test_flux_potential_trace_on_tips_and_internal_boundaries(base_discr: str):
     trace is equal to the pressure in the adjacent cell.
 
     """
-    model = DiffTpfaFractureTipsInternalBoundaries({"base_discr": base_discr})
+    model = DiffTpfaFractureTipsInternalBoundaries(
+        {"base_discr": base_discr, "times_to_export": []}
+    )
     model.prepare_simulation()
 
     mdg = model.mdg
diff --git a/tests/numerics/nonlinear/test_line_search.py b/tests/numerics/nonlinear/test_line_search.py
index c38500202d..0a677cb492 100644
--- a/tests/numerics/nonlinear/test_line_search.py
+++ b/tests/numerics/nonlinear/test_line_search.py
@@ -70,7 +70,7 @@ def test_line_search():
     checks are made on the results or that the line search does what it is supposed to.
 
     """
-    model = ConstraintFunctionsMomentumBalance()
+    model = ConstraintFunctionsMomentumBalance({"times_to_export": []})
     solver_params = {
         "nonlinear_solver": ConstraintLineSearchNonlinearSolver,
         "Global_line_search": True,
diff --git a/tests/numerics/nonlinear/test_nonlinear_solvers.py b/tests/numerics/nonlinear/test_nonlinear_solvers.py
index dc9652f195..44969e8977 100644
--- a/tests/numerics/nonlinear/test_nonlinear_solvers.py
+++ b/tests/numerics/nonlinear/test_nonlinear_solvers.py
@@ -40,7 +40,7 @@ def test_nonlinear_iteration_count():
     iteration count matches the pre set value after convergence is obtained.
 
     """
-    model = NonlinearSinglePhaseFlow()
+    model = NonlinearSinglePhaseFlow({"times_to_export": []})
     model.expected_number_of_iterations = 3
     pp.run_time_dependent_model(model, {})