diff --git a/integration-tests/model/test_prescribed_transport.py b/integration-tests/model/test_prescribed_transport.py
index cf31b2bfc..30cbc0292 100644
--- a/integration-tests/model/test_prescribed_transport.py
+++ b/integration-tests/model/test_prescribed_transport.py
@@ -1,11 +1,13 @@
 """
-This tests the prescribed wind feature of the PrescribedTransport timestepper.
+This tests the prescribed wind feature of the PrescribedTransport and
+SplitPresribedTransport (with no physics schemes) timesteppers.
 A tracer is transported with a time-varying wind and the computed solution is
 compared with a true one to check that the transport is working correctly.
 """
 
 from gusto import *
 from firedrake import sin, cos, norm, pi, as_vector
+import pytest
 
 
 def run(timestepper, tmax, f_end):
@@ -13,7 +15,8 @@ def run(timestepper, tmax, f_end):
     return norm(timestepper.fields("f") - f_end) / norm(f_end)
 
 
-def test_prescribed_transport_setup(tmpdir, tracer_setup):
+@pytest.mark.parametrize('timestep_method', ['prescribed', 'split_prescribed'])
+def test_prescribed_transport_setup(tmpdir, tracer_setup, timestep_method):
 
     # Make domain using routine from conftest
     geometry = "slice"
@@ -31,11 +34,19 @@ def u_evaluation(t):
                           sin(2*pi*t/setup.tmax)*sin(pi*z)])
 
     transport_scheme = SSPRK3(domain)
-    transport_method = DGUpwind(eqn, 'f')
 
-    timestepper = PrescribedTransport(eqn, transport_scheme, setup.io,
-                                      transport_method,
-                                      prescribed_transporting_velocity=u_evaluation)
+    if timestep_method == 'prescribed':
+        transport_method = DGUpwind(eqn, 'f')
+        timestepper = PrescribedTransport(eqn, transport_scheme, setup.io,
+                                          transport_method,
+                                          prescribed_transporting_velocity=u_evaluation)
+    elif timestep_method == 'split_prescribed':
+        transport_method = [DGUpwind(eqn, 'f')]
+        timestepper = SplitPrescribedTransport(eqn, transport_scheme, setup.io,
+                                               transport_method,
+                                               prescribed_transporting_velocity=u_evaluation)
+    else:
+        raise NotImplementedError
 
     # Initial conditions
     timestepper.fields("f").interpolate(setup.f_init)
diff --git a/integration-tests/physics/test_precipitation.py b/integration-tests/physics/test_precipitation.py
index 4c4a6441b..6f1b340e7 100644
--- a/integration-tests/physics/test_precipitation.py
+++ b/integration-tests/physics/test_precipitation.py
@@ -47,12 +47,12 @@ def setup_fallout(dirname):
 
     # Physics schemes
     rainfall_method = DGUpwind(eqn, 'rain', outflow=True)
-    physics_parametrisations = [Fallout(eqn, 'rain', domain, rainfall_method)]
+    physics_parametrisations = [(Fallout(eqn, 'rain', domain, rainfall_method), SSPRK3(domain))]
 
     # build time stepper
     scheme = SSPRK3(domain)
-    stepper = PrescribedTransport(eqn, scheme, io, transport_method,
-                                  physics_parametrisations=physics_parametrisations)
+    stepper = SplitPrescribedTransport(eqn, scheme, io, transport_method,
+                                       physics_schemes=physics_parametrisations)
 
     # ------------------------------------------------------------------------ #
     # Initial conditions
diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py
new file mode 100644
index 000000000..f49f87b62
--- /dev/null
+++ b/integration-tests/physics/test_terminator_toy.py
@@ -0,0 +1,110 @@
+"""
+This tests the terminator toy physics scheme that models the interaction
+of two chemical species through coupled ODEs.
+"""
+
+from gusto import *
+from firedrake import IcosahedralSphereMesh, Constant, cos, \
+    sin, SpatialCoordinate, Function, max_value, as_vector, \
+    errornorm, norm
+import numpy as np
+
+
+def run_terminator_toy(dirname):
+
+    # ------------------------------------------------------------------------ #
+    # Set up model objects
+    # ------------------------------------------------------------------------ #
+
+    # A much larger timestep than in proper simulations, as this
+    # tests moving towards a steady state with no flow.
+    dt = 50000.
+
+    # Make the mesh and domain
+    R = 6371220.
+    mesh = IcosahedralSphereMesh(radius=R,
+                                 refinement_level=2, degree=2)
+
+    # get lat lon coordinates
+    x = SpatialCoordinate(mesh)
+    lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2])
+
+    domain = Domain(mesh, dt, 'BDM', 1)
+
+    # Define the interacting species
+    X = ActiveTracer(name='X', space='DG',
+                     variable_type=TracerVariableType.mixing_ratio,
+                     transport_eqn=TransportEquationType.advective)
+
+    X2 = ActiveTracer(name='X2', space='DG',
+                      variable_type=TracerVariableType.mixing_ratio,
+                      transport_eqn=TransportEquationType.advective)
+
+    tracers = [X, X2]
+
+    # Equation
+    V = domain.spaces("HDiv")
+    eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V)
+
+    output = OutputParameters(dirname=dirname+"/terminator_toy",
+                              dumpfreq=10)
+    io = IO(domain, output)
+
+    # Define the reaction rates:
+    theta_c = np.pi/9.
+    lamda_c = -np.pi/3.
+
+    k1 = max_value(0, sin(theta)*sin(theta_c) + cos(theta)*cos(theta_c)*cos(lamda-lamda_c))
+    k2 = 1
+
+    physics_schemes = [(TerminatorToy(eqn, k1=k1, k2=k2, species1_name='X',
+                        species2_name='X2'), BackwardEuler(domain))]
+
+    # Set up a non-divergent, time-varying, velocity field
+    def u_t(t):
+        return as_vector([Constant(0)*lamda, Constant(0)*lamda, Constant(0)*lamda])
+
+    X_T_0 = 4e-6
+    X_0 = X_T_0 + 0*lamda
+    X2_0 = 0*lamda
+
+    transport_scheme = SSPRK3(domain)
+    transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')]
+
+    stepper = SplitPrescribedTransport(eqn, transport_scheme, io,
+                                       spatial_methods=transport_method,
+                                       physics_schemes=physics_schemes,
+                                       prescribed_transporting_velocity=u_t)
+
+    stepper.fields("X").interpolate(X_0)
+    stepper.fields("X2").interpolate(X2_0)
+
+    stepper.run(t=0, tmax=10*dt)
+
+    # Compute the steady state solution to compare to
+    steady_space = domain.spaces('DG')
+    X_steady = Function(steady_space)
+    X2_steady = Function(steady_space)
+
+    r = k1/(4*k2)
+    D_val = sqrt(r**2 + 2*X_T_0*r)
+
+    X_steady.interpolate(D_val - r)
+    X2_steady.interpolate(0.5*(X_T_0 - D_val + r))
+
+    return stepper, X_steady, X2_steady
+
+
+def test_terminator_toy_setup(tmpdir):
+    dirname = str(tmpdir)
+    stepper, X_steady, X2_steady = run_terminator_toy(dirname)
+    X_field = stepper.fields("X")
+    X2_field = stepper.fields("X2")
+
+    print(errornorm(X_field, X_steady)/norm(X_steady))
+    print(errornorm(X2_field, X2_steady)/norm(X2_steady))
+
+    # Assert that the physics scheme has sufficiently moved
+    # the species fields near their steady state solutions
+    assert errornorm(X_field, X_steady)/norm(X_steady) < 0.4, "The X field is not sufficiently close to the steady state profile"
+    assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.4, "The X2 field is not sufficiently close to the steady state profile"
diff --git a/integration-tests/physics/test_wind_drag.py b/integration-tests/physics/test_wind_drag.py
index 6bf107a9d..c6aec9bc4 100644
--- a/integration-tests/physics/test_wind_drag.py
+++ b/integration-tests/physics/test_wind_drag.py
@@ -37,7 +37,7 @@ def run_wind_drag(dirname, implicit_formulation):
     eqn = CompressibleEulerEquations(domain, parameters)
 
     # I/O
-    output = OutputParameters(dirname=dirname+"/surface_fluxes",
+    output = OutputParameters(dirname=dirname+"/wind_drag",
                               dumpfreq=1,
                               dumplist=['u'])
     io = IO(domain, output)