diff --git a/integration-tests/physics/test_terminator_toy.py b/integration-tests/physics/test_terminator_toy.py
new file mode 100644
index 000000000..39e85d938
--- /dev/null
+++ b/integration-tests/physics/test_terminator_toy.py
@@ -0,0 +1,111 @@
+"""
+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, ge, le, exp, cos, \
+    sin, conditional, interpolate, SpatialCoordinate, VectorFunctionSpace, \
+    Function, assemble, dx, FunctionSpace, pi, max_value, acos, 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 = 100000.
+
+    # Make the mesh and domain
+    R = 6371220.
+    mesh = IcosahedralSphereMesh(radius=R,
+                                 refinement_level=3, degree=2)
+    
+    domain = Domain(mesh, dt, 'BDM', 1)
+    
+    # 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))
+      
+    transport_scheme = SSPRK3(domain)
+    transport_method = [DGUpwind(eqn, 'X'), DGUpwind(eqn, 'X2')]
+    
+    X_0 = 4e-6 + 0*lamda
+    X2_0 = 0*lamda
+    
+    transport_scheme = SSPRK3(domain)#, limiter=MixedLimiter)
+    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)
+    
+    #True fields to compare to
+    
+    stepper.run(t=0, tmax=10*dt)
+    
+    steady_space = domain.spaces('DG')
+    X_steady = Function(steady_space)
+    X2_steady = Function(steady_space)
+    
+    X_T_0 = 4e-6
+    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")
+    
+    #Compute 
+    assert errornorm(X_field, X_steady)/norm(X_steady) < 0.25, "The X field is not sufficiently close to the steady state profile"
+    assert errornorm(X2_field, X2_steady)/norm(X2_steady) < 0.25, "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)